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Abstract 

This  thesis  describes  the  development  of  a  series  of  models  utilizing  the 
commercial  finite  element  suite  ABAQUS  specifically  to  apply  towards  the  study  of 
biological  tissue.  The  end  goal  is  to  be  able  to  obtain  the  material  properties  of  the 
Manducca  Sexta,  a  biological  inspiration  for  flapping  wing  micro-air  vehicles. 

Two  finite  element  models  were  used  to  analyze  the  results  of  two  prior  studies  of 
other  researchers.  A  flat  punch  elastic  model  examined  boundary  effects  and  confirmed 
that  the  point  of  indentation  was  far  enough  removed  from  the  boundary.  The 
hyperelastic  spherical  indentation  experiment  examined  the  effects  of  coefficient  of 
friction  on  the  indentation.  Another  algorithm  was  reproduced  to  analyze  the  elastic, 
power  law-hardening  properties  of  a  wide  range  of  material  properties. 

A  nanoindentation  system  was  used  to  investigate  the  modulus  of  the  M.  Sexta. 
Due  to  instrument  limitations,  useful  data  was  not  able  to  be  collected.  An  upper  bound 
on  the  modulus  was  analytically  established  at  approximately  3  MPa  using  the  noise  level 
of  the  equipment.  A  uniaxial  tension  test  of  the  M.  Sexta  was  used  to  obtain  a  reported 
initial  modulus  of  elasticity  values  of  343  kPa. 
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APPLICATION  OF  FINITE  ELEMENT  TO  EVALUATE  MATERIAL  WITH 
SMALL  MODULUS  OF  ELASTICITY 

I.  Introduction 


1.1  Objective 

This  study  describe  a  series  of  models  utilizing  the  commercial  finite  element 
suite  ABAQUS  which  allow  hyperelastic  materials  to  be  considered  for  an  application 
towards  soft  biological  tissue.  Specifically,  this  present  study  attempts  to  apply  the 
models  developed  within  this  study  to  evaluate  the  material  properties  of  the  muscles  of 
the  Manducca  Sexta.  This  insect,  known  as  the  hawkmoth,  is  of  interest  for  study  as  a 
biological  inspiration  for  flapping  wing  micro-air  vehicles  [1].  Attempts  to  model  the 
entire  moth  have  been  made  by  Demasi  [2];  however,  deficiencies  remain  in  the  accuracy 
of  the  input  material  to  that  model.  The  present  study  attempted  to  correct  that  deficiency 
through  finite  element  modeling  and  experimentation  with  nanoindentation,  along  with  a 
uniaxial  tensile  test.  As  part  of  the  study,  the  structure  of  the  muscle  and  information  on 
the  mechanics  of  insect  flight  are  discussed.  Also,  to  ensure  full  understanding  of  the 
process  on  nanoindentation,  models  using  standard  engineering  materials  such  as 
aluminum  are  examined  as  these  materials  have  better-characterized  properties  than 
hyperelastic  materials. 

1.2  Research  focus 

Nanoindentation  is  a  technique  that  can  be  used  to  analyze  problems  on  a  small 
length  scale  problems.  It  is  a  popular  method  for  evaluating  the  mechanical  properties  of 
materials  and  structures,  including  elastic  modulus,  yield  strength,  hardening  coefficient, 
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residual  stress,  fracture  toughness,  and  viscoelastic  behavior.  In  addition  to  its  usefulness 
in  examining  problems  on  a  small  length  scale,  the  preparation  required  for  these  tests 
can  be  less  intensive  than  traditional  tests,  such  as  uniaxial  or  biaxial  tests. 

Typically,  commercial  indenters  are  set  up  as  “black-box”  instruments  that 
automatically  calculate  the  elastic-plastic  properties  of  the  material.  However,  many 
materials  do  not  conform  to  the  elastic-plastic  theories  and  additional  analysis  of  the 
force-indentation  data  is  required.  Finite  element  (FE)  modeling  can  aid  in  this  additional 
analysis  by  correlating  the  shape  of  the  indentation  curve  with  the  desired  properties. 

The  hawkmoth  muscles,  called  the  dorsal  longitudinal  muscles  (DLMs)  and 
dorsal  ventral  muscles  (DVMs),  can  be  described  as  behaving  as  an  orthotropic, 
viscoelastic  material.  Spherical  indentation  experiments  (utilizing  the  developed  finite 
element  model)  were  attempted  to  examine  the  transverse  elastic  properties  and  the  finite 
element  models  were  use  to  show  why  they  failed.  Uniaxial  tension  tests  were  used  to 
attempt  to  obtain  the  longitudinal  properties  of  the  specimen. 

1.3  Methodology 

All  finite  element  modeling  was  accomplished  using  the  ABAQUS  finite  element 
program  from  Dassault  Systems.  FE  models  were  developed  then  validated  with 
experimental  data  and  closed  form  analytic  solutions,  when  available.  Experiments 
involving  nanoindentation  were  attempted  using  the  MTS  Nanoindenter  XP.  This 
equipment  was  used  with  the  permission  of  the  Air  Force  Research  Laboratories  (AFRL) 
Materials  and  Manufacturing  Directorate.  All  uniaxial  tensile  tests  were  performed  using 
the  MTS  Nano  Universal  Testing  Machine.  Moth  specimens  used  to  obtain  the  muscle 
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samples  were  provided  by  Case  Western  University.  More  information  on  the  equipment 
used  in  this  study  and  instructions  for  raising  the  moths  can  be  found  in  the  Appendix. 


1.4  Literature  Review 

1.4.1  Indentation  Testing  Introduction 

In  this  section,  a  brief  overview  of  nanoindentation  is  given  to  introduce  what  it  is 
and  what  a  typical  experiment  involves. 

In  nanoindentation  testing,  a  probe  is  pressed  into  a  material  surface  under  either 
load-  or  displacement-control.  The  load  and  displacement  data  can  be  used  to  determine 
various  material  properties  of  the  sample.  Useful  results  from  nanoindentation 
experiments  require  measurement  of  extremely  small  forces  and  displacements  with  great 
accuracy  and  sensitivity.  A  simple  indentation  diagram  (a)  and  a  typical  load- 
displacement  curve  (b)  can  be  seen  in  Figure  1 . 


Figure  1:  (a)  Simple  Diagram  of  Indentation  (b)  Example  Load-Displacement  Plot 
for  a  Standard  Instumented  Indentation 
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Oliver  and  Pharr  pioneered  the  effective  use  of  a  small-scale  version  of  the 
standard  indentation  test  for  an  elastic-plastic  material  [3],  Nanoindentation  had  been 
around  for  a  decade  at  that  point,  but  the  results  had  not  been  as  precise  or  repeatable  as 
other,  more  established  methods  for  testing  material  properties.  One  of  the  most 
significant  contributions  of  their  work  was  the  realization  that  the  unloading  curve  of  the 
load-displacement  plot  was  not  linear,  as  had  previously  been  assumed.  The 
mathematical  relationships  used  to  determine  the  properties  are  discussed  in  depth  in 
Chapter  2. 

As  mentioned,  two  types  of  indenters  were  examined  in  this  study:  cylindrical  and 
spherical  (Figure  2).  The  geometry  of  the  probe  is  important  because  the  shape 
determines  the  deformation  profiles  obtained  during  the  test.  The  cylindrical,  flat-tips  are 
convenient  because  they  have  a  constant  contact  area  with  the  sample.  This  simplifies 
much  of  the  analysis,  however,  for  small-scale  tests  it  can  be  difficult  to  align  the  probe 
surface  with  the  surface  of  the  sample.  Spherical  tips  are  advantageous  in  that  they  delay 
the  immediate  onset  of  plastic  deformation.  This  is  one  of  the  reasons  this  indenter  tip  is 
popular  with  evaluation  of  biological  tissues. 
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(a)  (b) 

Figure  2:  Cartoon  Schematic  of  Common  Indenter  Geometries,  (a)  Cylindrical  (b) 

Spherical,  Indenting  a  Flat  Half-space 


1.4.2  Elastic-Plastic  Spherical  Finite  Element  Model 

A  single  sharp  indentation  cannot  yield  a  unique  solution  [4]  through  the  standard 
indentation  technique.  For  a  power-law  hardening  material  with  properties  ay  (yield 
strength),  E,  v,  and  n  (work-hardening  coefficient),  only  2  of  the  properties  can  be 
determined  while  the  other  2  have  to  be  known  a  priori.  However,  researchers  have  been 
able  to  take  advantage  of  a  concept  called  representative  strain  to  simplify  the  problem 
using  finite  element  analysis.  [5]  Zhao,  et  al.  used  spherical  indentation  along  with 
representative  strain  to  determine  the  unique  properties  from  one  indentation.  From  that 
one  indentation  they  selected  two  points  on  the  force-displacement  loading  curve  and  the 
contact  stiffness  (slope  of  the  unloading  curve)  to  produce  3  equations  to  be  able  to  solve 
for  CTy,  E,  and  n  (v  was  found  to  vary  little  for  the  analysis  and  was  assumed  to  be  0.3). 
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p 


Figure  3:  Representative  Location  of  2  Points  Required  for  the  Elastic-Plastic 
Spherical  Model  Shown  on  a  Generic  Loading  and  Unloading  Diagram 

The  goal  of  this  study  was  to  utilize  the  algorithm  developed  by  Zhao  for  a  wide  range  of 
materials.  Possible  applications  to  the  exoskeleton  of  the  hawkmoth  are  examined.  The 
formulation  of  these  equations  is  discussed  in  Chapter  2  and  the  finite  element  model  is 
reproduced  in  Chapter  3. 

1.4.3  Biomechanics  Introduction 

Experiments  to  determine  the  material  properties  of  muscles  and  other  soft  tissues 
have  been  documented  for  centuries.  Many  of  the  famous  physical  relationships  that  are 
used  in  engineering  today  had  their  origins  in  scientists  investigating  biomechanics. 
Leonhard  Euler  was  the  first  to  examine  propagation  of  pulse  waves  in  arteries.  Thomas 
Young  (of  Young’s  Modulus  fame)  studied  the  formation  of  human  voice  and  connected 
it  with  the  elasticity  of  the  vocal  cords.  Poiseulle  studied  the  pressure  in  the  aortas  of 
dogs  which  led  to  the  establishment  of  the  no-slip  condition  in  pipe  flow.  [6] 
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The  determination  of  mechanical  properties  of  soft  tissues  in  biomechanics  is 
often  challenging.  Isolating  the  tissue  from  the  subject  for  testing  is  often  very  difficult. 
The  small  size  of  the  samples  and  the  need  to  maintain  hydration  and  in  normal  living 
conditions  can  cause  complications  as  well.  In  addition,  biological  tissues  are  often 
nonlinear  and  dependent  on  the  history  of  the  loading  and  unloading  cycle.  More  of  these 
characteristics  are  discussed  in  Chapter  2. 

Uniaxial,  along  with  biaxial  (longitudinal  and  transverse  direction),  tension  tests 
are  some  of  the  more  common  methods  to  characterize  the  properties  of  soft  tissues.  The 
tension  test  is  often  used  as  baseline  from  which  other  experiments  branch  off.  Early 
experiments  by  Wertheim  (1847)  showed  the  elastic  nature  of  soft  tissue.  Hill  (1938) 
applied  data  from  uniaxial  tests  for  his  first  models  of  the  contraction  of  muscles  [5], 
Other  tests  of  note:  Moss  and  Halpem  (1977)  [7]  determined  the  viscous  and  elastic 
properties  of  resting  frog  muscle;  Van  Locke  et  al.,  (2006)  [8]  examined  the  compressive 
behavior  of  muscle;  and  Lally  et  al.,  (2004)  [9]  studied  the  effects  of  biaxial  and  uniaxial 
tension  on  pig  artery 

The  popularity  of  indentation  tests  for  biological  materials  has  grown  quickly  in 
the  last  decade  as  the  machinery  and  techniques  required  for  the  tests  has  advanced.  In 
the  following  few  sections,  several  different  indentation  and  uniaxial  experiments  with 
soft  tissues  are  reviewed.  The  techniques  used  in  these  examples  were  the  basis  for  the 
finite  element  models  and  the  attempted  experiments  in  this  present  study  are  described 
in  more  detail  in  Chapters  2  and  3. 
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1.4.4  Soft  Tissue  Indentation  FE  and  Experimentation 

Conducting  a  nanoindentation  experiment  with  a  biological  tissue  sample  is  a 
very  challenging  task.  There  are  many  factors  that  can  influence  the  experiment.  The 
factors  influencing  sample  preparation  and  tip  preparation  described  in  this  section  are 
selected  from  the  Handbook  of  Nanoindentation  with  Biological  Applications  [10],  [11], 

One  of  the  most  influential  factors  in  indenting  biological  tissue  is  maintaining  its 
hydration.  Soft  tissues  are  made  up  of  mainly  water  and  exposure  to  air  results  in  tissue 
desiccation.  This  desiccation  can  change  the  material  properties  of  the  sample. 

Biological  tissue  samples  also  have  a  surface  roughness  that  would  otherwise  be 
removed  through  a  polishing  process  for  indention  of  a  standard  engineering  material. 
This  technique  is  not  available  when  testing  soft  tissues  as  their  microstructures  could  be 
substantially  disrupted.  This  roughness  influences  the  tip  selection.  The  low  modulus  of 
the  tissue  normally  requires  the  use  of  cylindrical  flat  punch  or  spherical  tips.  However, 
the  spherical  tip  allows  for  some  inaccuracy  in  the  approach  to  the  sample  and  is 
therefore  often  used  for  testing  soft  tissues  with  irregular  surfaces. 

Also,  many  soft  tissues  are  not  isotropic.  The  anisotropic  nature  doesn’t  meet 
many  of  the  assumptions  for  standard  indentation  theory.  Since  the  muscle  sample  has 
properties  that  are  transversely  orthotropic,  the  results  of  this  test  would  produce  an 
indentation  elastic  modulus  that  is  a  function  of  the  transverse  and  longitudinal  moduli, 
weighted  in  the  direction  of  the  indentation,  the  transverse  direction.  The  materials  also 
have  large  displacements  for  a  given  load  due  to  their  hyperelastic  nature. 

Instrumentation  also  presents  limitations  on  testing  soft  tissues.  Most  commercial 
instrumented  nanoindentation  devices  are  calibrated  for  testing  materials  with  a  modulus 
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in  the  GPa  or  MPa  range.  Most  biological  soft  materials  have  elastic  moduli  in  the  range 
of  tens  to  hundreds  of  kPa’s  and  can  present  problems  with  finding  the  surface  of  the 
material.  As  is  discussed  in  Chapter  5,  the  resolution  of  the  chosen  instrumentation 
proved  to  be  the  source  of  the  difficulty  in  completing  this  experiment. 

To  aid  in  the  development  of  the  procedures  for  the  experiment,  several 
previously-conducted  experiments  and  their  associated  finite  element  models  (if  used) 
were  examined.  The  first  indentation  test  described  in  this  section  comes  from  a  group  at 
the  National  Institute  of  Health  who  investigated  the  spherical  indentation  of  soft  matter 
in  the  hyperelastic  regime.  [12]  Researchers  conducting  spherical  nanoindentation 
experiments  for  Hookean  elastic  materials  compare  their  results  to  the  classical  Hertz 
solution  for  a  sphere  impacting  a  flat  plane.  This  study  focuses  on  developing  a 
relationship  for  hyperelastic  materials  through  finite  element  modeling  with  several 
different  hyperelastic  strain  energy  potential  functions.  Lin,  et  al.  compared  the  resulting 
functions  to  synthetic  gels  and  mouse  cartilage.  The  mouse  cartilage  test  yielded  a  shear 
modulus  of  14.3  kPa  (p)  and  a  fitting  parameter  (a)  of  7.3.  The  derivations  of  Hertzian 
contact  and  the  hyperelastic  relations  are  discussed  in  Chapter  2.  The  hyperelastic  finite 
element  model  used  in  this  study  was  based  on  the  Lin  paper. 

Another  indentation  test  that  was  examined  was  one  by  V.T.  Nayar  et  al.  [13]  and 
examined  porcine  sclera  (the  white  of  the  pig’s  eye).  The  results  from  this  study  are  not 
useful  for  comparison  to  the  muscles  of  the  hawkmoth  as  the  sclera  has  a  planar  isotropic 
structure.  However,  the  methods  used  in  this  study  were  helpful.  The  sclera 
(approximately  1-1.2  mm  thick)  was  removed  from  the  pig  in  approximately  1  cm 
squares  and  secured  to  a  glass  slide  with  cyanoacrylate  (superglue).  A  shallow  well  was 
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built  up  around  the  outside  of  the  sample  to  form  a  ring.  This  ring  (shown  in  Figure  4) 
was  then  filled  with  saline  solution  to  maintain  hydration  of  the  sample.  Testing  was 
accomplished  with  an  80  pm  cylindrical  punch.  Load  controlled  indentations  were 
conducted  to  375,  750,  and  1500  pN.  Shear  modulus  reported  was  approximately  30 
kPa. 


Figure  4:  Experimental  Setup  for  Porcine  Sclera  Indentation  Experiment  From 

Nayer  et  al.  [13] 


A  third  indentation  test  was  accomplished  with  the  skinned  cardiac  muscle  fibers 
of  a  cow.  The  cardiac  muscles  were  isolated  from  the  adult  cow,  rinsed  thoroughly,  and 
secured  in  a  mica  sheet.  This  experiment  used  atomic  force  microscopy  (AFM)  so  the 
experimental  protocols  would  not  be  comparable.  Their  results  from  this  technique  were 
on  the  order  of  20  kPa  for  the  shear  modulus.  [14] 
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1.4.5  Soft  Tissue  Uniaxial  Experimentation 

The  first  uniaxial  tension  study  discussed  in  depth  here  from  Laure-Lise  Gras,  et 
al.  [16],  examined  the  passive  response  of  human  stemocleidomastoideus  muscle  (located 
in  the  neck)  to  tension.  Their  study  was  performed  in  vitro  therefore  the  muscle  samples 
were  removed  from  cadavers  and  placed  into  a  uniaxial  testing  machine.  Boundary 
conditions  were  maintained  by  removing  portions  of  the  jaw  and  allowing  the  entire 
muscle  to  be  tested.  Typical  specimens  were  134  mm  long  with  a  cross-sectional  area  of 
300  mm  .  3D  reconstructions  can  be  seen  in  Figure  5. 

A  B  C 

II 

Figure  5:  Muscle  3D-reconstruction  and  mesh.  (A)  Example  of  a  3D  reconstruction. 

(B)  Its  associated  finite-element  model.  (C)  Superposition  of  the  3D  reconstruction 

and  the  finite-  element  model. 

In  order  to  prevent  the  specimen  from  desiccating,  Laure-Lise  Gras,  et  al.  would 
moisten  the  surface  of  the  muscles  regularly  with  a  saline  solution.  After 
preconditioning,  the  specimens  were  subjected  to  a  maximum  deformation  of  15%  at  a 
strain  rate  of  0.00125s'1.  After  assuming  an  incompressible  (v  =  0.5)  and  isochoric 
material,  results  were  fitted  to  an  exponential  hyperelastic  form  from  Stem-Knudsen  and 
also  to  the  Ogden  hyperelastic  constitutive  law.  These  results  were  compared  to  finite 
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element  models  obtained  by  iterating  the  parameters  of  the  function.  For  their  results, 
they  were  able  to  obtain  a  value  for  shear  modulus  of  37  kPa. 

Another  uniaxial  tension  test  was  accomplished  by  Calvo,  et  al.  [17]  only  their 
specimen  of  interest  was  rat  muscle  (Figure  6).  The  tendons  of  the  rat  were  also  examined 
my  Calvo,  et  al.  but  this  is  not  discussed  here.  The  samples  studied  by  Calvo,  et  al.  were 
much  closer  to  the  size  of  the  muscle  samples  of  the  hawkmoth.  Specimens  averaged  a 
length  of  6.6  mm  and  a  cross  sectional  area  of  1 .8  mm2.  The  samples  were  glued  to 
pieces  of  sandpaper  to  be  placed  into  the  grips  of  a  displacement  controlled  microtester. 

A  testing  velocity  of  0.2*L/100  mm  min'1  was  used,  which  for  the  average  muscle  length 
corresponded  to  3.3e-5  s-1  strain  rate,  which  is  slower  than  the  previous  study.  Sample 
hydration  was  maintained  through  a  cooled  ultrasonic  humidifier.  The  Calvo  study  used 
a  modified  form  of  the  Weiss  strain-energy  density  function  to  fit  the  experimental  data. 
Again,  isochoric  and  incompressible  assumptions  were  made.  One  thing  that  stands  out 
in  their  test  is  that  they  make  no  mention  of  preconditioning  the  sample  prior  to  testing. 
This  differs  from  standard  practice  outlined  and  could  be  a  source  of  variability  in  their 
data  due  to  viscoelastic  effects.  [5] 


Figure  6:  Rat  Muscle  Sample  with  Sandpaper  Grips  from  Calvo,  et  al. 
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A  third  uniaxial  tensile  actually  uses  M.  Sexta  as  a  subject,  only  with  a  very 
important  distinction.  The  group  from  Tufts  University  [18]  studied  the  muscle  of  the 
species  when  it  was  in  its  caterpillar  state  prior  to  its  metamorphosis.  An  initial  thought 
would  be  that  a  direct  comparison  could  be  made  between  this  muscle  and  the  muscle  of 
the  flying  insect.  However,  these  muscles  are  used  for  entirely  different  tasks.  The 
muscle  of  the  caterpillar  undergoes  strains  of  30%  in  approximately  1  second.  This 
differs  greatly  from  that  on  the  moth  where  strains  of  7%  in  around  0.018  s  are  the  norm. 
The  techniques  used  in  the  Tufts  study,  however,  are  useful  for  application  to  the  current 
study.  Both  passive  and  active  muscles  were  studied  which  provides  a  useful  comparison 
between  the  two  states. 

Muscle  samples  from  the  caterpillars  in  the  Tufts  study  were  approximately  4-5 
mm  long.  Cross  sectional  area  was  not  reported.  Muscles  were  pinned  by  the  attached 
cuticle  at  each  end  in  a  horizontal  bath  of  saline  to  prevent  dehydration.  One  end  of  the 
muscle  was  secured  by  a  hook  to  a  displacement  controlled  testing  machine.  The 
samples  were  preconditioned  for  between  6-10  cycles  to  remove  the  hysteresis.  The 
Dorfman  model  applied  a  modified  pseudo-elastic  model  from  Dorfmann  and  Ogden 
(2003).  Shear  modulus  reported  in  the  0.78  kPa  range  was  much  lower  than  the  other 
experiments. 

1.4.6  Biological  Specimen 

In  order  to  understand  the  reason  for  investigating  the  hawkmoth  muscle,  it  is 
useful  to  know  more  about  the  species  and  how  the  muscles  induce  flight.  This  is 
described  in  this  section. 
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The  biological  species  chosen  to  be  studied  for  the  experimentation  portion  of  this 
study  is  the  Manduca  Sexta.  The  species,  which  will  be  referred  to  as  the  hawkmoth  for 
the  remainder  of  the  study,  is  shown  in  Figure  7.  The  size,  weight,  and  flight 
performance  characteristics  of  the  hawkmoth  made  it  an  ideal  candidate  of  study  as  a 
natural  MAV.  The  hawkmoth  weighs  approximately  1 .5  to  2  g  and  has  a  wingspan 
between  9.5  and  12  cm.  [19]  As  such,  it  is  one  of  the  larger  flying  insects  in  nature. 


Figure  7:  A  Natural  Flapping-wing  MAV,  Manduca  Sexta 


The  hawkmoth  anatomy  can  be  divided  into  4  main  parts:  head,  wings,  thorax  and 
abdomen.  The  head  contains  the  primary  nervous  system  control  organ,  two  eyes  and 
antennas,  and  a  coiled  proboscis  for  feeding  on  plant  nectar.  The  abdomen  contains 
many  of  the  body’s  organs  for  digestion  and  reproduction.  [20]  Two  sets  of  wings,  the 
forewings  (larger  and  towards  head)  and  hind  wings  (smaller  and  closer  to  abdomen), 
consist  of  a  thin,  flexible  membrane  overlaying  a  network  of  rigid  veins.  The  wing’s 
membrane  is  covered  with  scales  that  are  used  for  camouflage  and  possibly  influence 
flow  patterns  during  flapping  [20].  The  thorax  is  located  at  the  intersection  of  the  other 
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main  parts  and  contains  the  mechanism  to  generate  and  control  flight.  A  cross-sectional 
diagram  of  the  thorax  can  be  seen  in  Figure  8.  Within  the  exoskeleton  of  the  thorax  are 
the  DVMs  and  DLMs  as  well  as  a  bundle  of  nerves  called  the  ganglia  that  act  as  a 
secondary  brain  to  control  the  movement  of  the  muscles.  On  top  of  the  thorax  is  a  much 
thicker  and  more  rigid  section  called  the  tergal  plate.  The  interaction  of  the  tergal  plate 
and  the  flight  muscles  results  in  the  wings  flapping  as  is  shown  in  the  next  section.  A 
simplified,  3-dimensional  diagram  of  the  thorax  and  the  flight  muscles  can  be  seen  in 
Figure  9. 


Figure  8:  Cross-sectional  diagram  of  thorax  highlighting  the  various  muscle  groups 
in  the  Hawkmoth  and  the  interaction  with  Exoskeleton 
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Figure  9:  3D  Diagram  of  Thorax  and  Flight  Muscles 

The  hawkmoth  is  known  for  its  impressive  ability  to  hover.  Its  rapid  wingbeat 
has  often  led  to  it  being  misidentified  as  a  hummingbird.  The  hawkmoth  accomplishes 
its  hover  through  what  is  known  as  synchronous,  indirect  flight.  During  synchronous 
flight,  for  every  pulse  from  a  neuron  there  is  a  contraction  of  the  muscles  and  one 
corresponding  flapping  motion  of  the  wings.  Asynchronous  flight  differs  from 
synchronous  flight  in  that  one  neuron  pulse  will  cause  multiple  contractions  of  the 
muscles  which  will  produce  multiple  cycles  of  flapping.  Synchronous  flight  is  common 
in  insects  with  flapping  frequency  below  100  Hz.  An  indirect  flapping  mechanism  is  one 
in  which  the  muscles  do  not  directly  attach  to  the  wings.  They  pull  on  the  exoskeleton 
which  then  translates  that  motion  into  flapping  via  a  hinge.  Direct  flight  is  produced  by 
the  flight  muscles  attaching  directly  to  the  wing. 

The  indirect  flight  of  the  hawkmoth  begins  with  the  DVMs.  The  DVMs  (shown 
in  Figure  10  in  the  brown)  are  oriented  vertically  and  at  an  angle.  These  muscles  contract 
and  compress  the  exoskeleton.  The  compression  causes  the  inner  portion  of  the  wing 
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hinge  to  drop  and  the  wing  to  rise.  To  produce  the  downstroke,  the  DLMs  (shown  in 
pink)  contract  and  cause  the  exoskeleton  to  bow  upward.  The  bowed  exoskeleton  raises 
the  hinge  and  lowers  the  wings.  The  DLMs  are  the  much  larger,  and  therefore  more 
powerful,  of  the  two  sets.  The  DLMs  relatively  powerful  downstroke  produces  the 
majority  of  the  lift.  It  has  been  shown  in  several  studies  of  the  moth  that  both  the 
upstroke  and  downstroke  will  produce  lift  via  changing  the  camber  of  the  wing.  The 
change  in  camber  is  passive;  the  moth  does  not  directly  control  it.  A  typical  hawkmoth 
will  flap  its  wings  at  a  frequency  around  20  or  25  Hz. 

■  DLM  ■  DVM 
#  Hinge 


Wins 


A 

Figure  10:  Diagram  of  Flight  Mechanics  of  Hawkmoth:  (A)  DVM 
Contraction/Upstroke;  (B)  DLM  Contraction/Downstroke  [23] 

The  most  basic  contractile  unit  of  the  muscle  is  the  sarcomere.  The  sarcomere  is 
a  region  of  interaction  of  two  myofilaments,  one  thick  and  one  thin.  The  thick  filaments 
are  myosin  molecules,  while  the  thin  ones  are  actin  molecules.  (Diagram  shown  in  Figure 
11)  The  sarcomere  is  approximately  2.5  pm  long,  with  the  exact  length  dependent  on  the 
force  acting  in  the  muscle  and  the  state  of  excitation.  [6]  When  the  muscle  receives  a 
signal  from  the  motoneuron,  the  molecules  attached  to  the  myosin  filaments  extend  out, 
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pulling  on  the  actin  filaments.  This  causes  the  two  filaments  to  slide  past  each  other  and 
contract  the  muscle.  Neither  the  myosin  or  actin  filament  shortens  during  the  contraction. 
[6] 


Actin 


Sarcomere 


t 

Myosin 


Figure  11:  Diagram  of  the  basic  functional  unit  of  Muscle  Tissue  (Sarcomere)  made 
up  of  Interlocking  Actin  and  Myosin  Filaments. 

The  rest  of  the  structural  arrangement  on  the  muscle  can  analogous  to  a  Russian 
nesting  doll.  The  myofilaments  are  bound  into  groups  of  called  muscle  fibrils.  Muscle 
fibers  contain  groups  of  these  fibrils.  All  of  the  muscle  fibers  that  are  innervated  from  a 
single  motoneuron  are  called  muscle  motor  units.  The  complete  muscle  (the  outermost 
doll)  is  made  up  of  all  the  motor  units.  As  mentioned  earlier,  in  the  hawkmoth  there  are 
two  sets  of  flight  muscles:  the  DVMs  and  DLMs.  Each  DLM  is  composed  of  5  motor 
units,  numbered  1  through  five  on  the  right  hand  side  of  Figure  12.  The  DVMs  contain  6 
motor  units  each;  however,  they  are  less  massive  than  the  powerful  DLMs. 
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Figure  12:  Side  View  Schematic  (Left)  and  CT  scan  (Right)  of  M.  Sexta  Primary 

Flight  Muscles 

Other  members  of  the  research  team  here  at  AFIT  have  analyzed  different  aspects 
of  this  species  physiology  to  investigate  these  flight  mechanics.  Most  recently  Major 
Ryan  O’Hara  researched  the  material  properties  and  structural  dynamics  of  the  forewing. 
[20]  Lt  Alex  Hollenbeck  [21]  and  Brian  Cranston  [22]  investigated  the  exoskeleton  of  the 
thorax  material  properties  and  how  compression  in  the  vertical  and  longitudinal 
directions  affected  the  flapping  motion  of  the  wings  and  related  it  to  power  output  from 
the  muscles.  Also,  Captain  Travis  Tubbs  explored  the  timing  of  the  muscle  neurons 
through  electromyography.  [23]  This  study  is  the  first  here  at  AFIT  to  attempt  to 
directly  look  at  the  material  properties  of  the  muscles. 

1.5  Research  Implications 

The  finite  element  models  developed  during  this  study  are  valuable  to  other 
researchers  here  at  AFIT.  The  finite  element  models  have  been  validated  using 
experimental  data  and  analytic  equations  so  there  will  be  considerable  time  savings  for 
another  student  and  aid  them  in  their  research. 
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One  possible  application  of  the  insight  gained  from  the  experimental  portion  of 
this  study  could  be  to  develop  artificial  muscles  that  could  power  an  artificial  flapping- 
wing  MAV.  Muscles  have  been  used  for  inspiration  in  design  since  the  1950s  and  60s. 
The  McKibben  pneumatic  artificial  muscles,  developed  by  the  Bridgestone  Rubber 
Company  of  Japan,  were  an  early  example.  [24]  These  devices  pumped  pressurized  air 
into  a  rubber  bladders  enclosed  in  a  mesh  shell  to  mimic  a  muscle  contraction.  These 
devices  produced  narrower  dynamic  range,  but  higher  force  intensity  than  natural 
muscles.  Applications  for  MAV  development  would  not  be  very  practical  due  to  the  size 
and  power  constraints  of  the  devices. 

More  recently,  new  electroactive  polymers  (EAPs)  have  been  developed  that  have 
the  properties  required  to  simulate  the  contraction  of  muscles  [25],  [26],  EAPs  are 
materials  which  can  change  shape  in  response  to  an  electrical  stimulus.  Since  the  1990s, 
the  strains  capable  of  being  produced  by  these  materials  have  increased  thereby 
expanding  their  usefulness.  Two  of  the  latest  advances  in  EAP  technology  within  the  last 
year  involve  carbon  nanotube  aerogels  [27]  and  telescopic  polymer  chains  [28], 

1.6  Thesis  Preview 

In  Chapter  2,  the  theory  behind  nanoindentation  is  outlined.  The  analytic 
solutions  for  indentations  into  an  elastic  material  are  outlined.  These  equations  are  used 
to  compare  against  the  finite  element  models  in  Chapter  3.  Additionally  the  constitutive 
equations  for  the  elastic,  hyperelastic,  and  elastic-plastic  material  properties  are  described 
as  well  as  the  material  structure  of  the  hawkmoth  muscle. 
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In  Chapter  3,  the  development  of  the  finite  element  models  is  described.  The 
boundary  conditions,  element  types,  and  the  model  type  are  discussed  for  each  of  the  2 
types  of  indenters  (spherical  and  cylindrical  flat  punch). 

In  Chapter,  4  the  hawkmoth  muscle  experimentation  methodology  is  discussed  for 
the  nanoindenter  and  the  uniaxial  tensile  experimentation.  The  process  of  dissection  for 
the  moths  is  also  outlined. 

In  Chapter  5,  the  3  models  are  used  in  scenarios  with  applications  to  testing 
biological  tissues.  Additionally,  problems  with  the  nanoindentation  experiment  are 
attempted  to  be  explained.  Lastly,  the  uniaxial  tension  test  results  is  reported  and 
discussed.  Chapter  6  summarizes  the  results  of  this  study  and  gives  recommendations  for 
future  research. 

II.  Theory 


2.1  Chapter  Overview 

The  purpose  of  this  Chapter  is  to  outline  the  analytic  solutions  of  the  interaction 
between  the  indenter  geometry  and  the  materials  in  the  models  developed.  These  theories 
are  the  foundation  for  the  models  of  nanoindentation  developed  in  Chapter  3.  The  first 
section  describes  a  fully  elastic  material  indented  by  a  probe  with  2  different  shapes 
(spherical  and  cylindrical  flat  punch).  In  the  second  section,  the  elastic  material  model 
material  properties  of  the  sample  are  changed  to  a  hyperelastic  material.  This  indenter 
model  is  compared  to  equations  derived  from  Lin  et  al  [12].  The  third  and  fourth 
sections  are  models  of  indentation  into  elastic-plastic  materials  by  a  sharp  and  spherical 
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indenter,  respectively.  The  spherical  indenter  model  was  developed  to  obtain  the  entire 
stress-strain  curve  from  a  single  indentation. 

Additionally,  as  mentioned  in  the  objective,  an  overview  of  the  muscle  material  of 
the  hawkmoth  is  laid  out.  That  final  section  outlines  the  models  of  the  tissue  structure. 
The  equations  used  in  a  uniaxial  tensile  test  are  also  described. 

2.2  FE  Models  Material  Theory 
2.2.1  Elastic  Indentation  Theory 

The  first  finite  element  model  described  in  this  research  study  is  a  linear  elastic, 

isotropic  material  that  behaves  according  to  the  following: 

su  1  -v  -v  0  0  0  (Tn 

£21  —v  1  —v  0  0  0  <j21 

£•33  1  -v  -v  1  0  0  0  C33 

2^3  ~  E  0  0  0  2(1  +  1/)  0  0  cr23 

2  £u  0  0  0  0  2(1  +  v)  0  cr13 

2  sl2_  0  0  0  0  0  2(l  +  v)__cr12 

where,  eq  and  o,j  are  the  strains  and  stresses  of  the  material,  and  E  and  v  are  the  elastic 
modulus  and  Poisson’s  ratio. 

The  analytic  solutions  to  be  used  to  compare  against  the  elastic  finite  element 
model  can  be  modeled  as  contact  between  two  elastic  bodies.  These  solutions,  which  can 
serve  as  a  reference  point  for  the  interaction,  were  first  studied  by  Hertz  in  1881. 
Equations  given  in  this  section  are  derived  from  Fischer-Cripps  Introduction  to  Contact 
Mechanics  [29],  For  his  formulations  Hertz  assumed: 
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i.  The  displacements  and  stresses  must  satisfy  the  differential 
equations  of  equilibrium  for  elastic  bodies  and  the  stresses 
must  vanish  at  a  great  distance  from  the  contact  surface. 

ii.  The  bodies  are  in  frictionless  contact. 

iii.  At  the  surface  of  the  bodies,  the  normal  pressure  is  zero 
outside  and  equal  and  opposite  inside  the  circle  of  contact. 

iv.  The  distance  between  the  surfaces  of  the  two  bodies  is  zero 
inside  and  greater  than  zero  outside  the  circle  of  contact. 

v.  The  integral  of  the  pressure  distribution  within  the  circle  of 
contact  with  respect  to  the  area  of  the  circle  of  contact  gives 
the  force  acting  between  the  two  bodies. 


The  elastic  modulus  of  the  contact  can  be  expressed  as  a  sum  of  the  two  bodies  by 
the  following  relation: 


(2.2) 


1  l-of  ,  l-^2 
ER  E,  es 


where  Er  is  the  reduced  modulus,  E,  and  Es  are  the  elastic  moduli  of  the  indenter  and 
sample,  respectively,  and  v;  and  vs  are  Poisson’s  ratios  of  the  indenter  and  sample, 
respectively.  If  the  indenter  is  assumed  to  have  a  much  larger  modulus  than  the  sample, 
the  equation  simplifies  to: 


(2.3) 


77  - _ 

R~i -t/ 


For  a  spherical  indenter  as  shown  in  Figure  13,  the  force-displacement  relationship  is: 
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(2.4) 


P  = 


R 


where  P  is  the  force  applied,  R  is  the  sphere  radius,  and  a  is  the  contact  radius.  For  a 

sphere,  the  contact  radius  can  be  expressed  as  a=yl~Rh ,  where  h  is  the  indentation  depth. 
Substituting  this  into  (2.4)  gives: 


(2.5) 


P  =  -E„R2h2 
3  * 


Figure  13:  Elastic  Spherical  Indentation 

In  addition  to  force-displacement  relationship,  Hertz  also  developed  equations  for  the 
stresses  from  the  indentation.  In  the  following  equation,  r  is  defined  as  a  point  along  the 
contact  between  the  two  bodies  that  originates  from  the  centerline  (i.e.  r=0  at  centerline). 


(2.6) 


-P  3  r2 

°'z=— tJi-— 

na  2  V  a 


For  the  spherical  indenter  geometry  presented  here,  the  indentation  stress  (or 
mean  pressure,  a*)  and  strain  (e*)  are  given  by 

.  P 

<7  =- 

(2.7) 


n  a 


£  =0.2 


a 


R 
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where  0.2  is  an  empirically  determined  constant  by  Tabor  (1951).  [12]  This  constant  has 
been  verified  by  other  investigations  as  well.  Substituting  the  force-displacement  into 

(2.7)  yields: 

(2.8)  a=—ERe 

3  n 

This  is  the  stress  strain  curve  for  a  linear  elastic  solid  indented  by  a  rigid  sphere. 

For  a  cylindrical,  flat  punch  as  shown  in  Figure  14,  the  force-displacement 
relationship  and  the  stress  field  equation  are: 


(2.9) 


P  —  2aERh 


(2.10) 


_  ~p  i  rz 

°”z  7ra2  2  V  a 2 


,  r  <a 


Figure  14:  Cylindrical  Flat  Punch  Indentation 


For  the  cylindrical  indenter  geometry  presented  here,  the  indentation  stress  (or 
mean  pressure,  a*)  and  strain  (s*)  are  given  by 


(2.11) 


* 

<J 


P 

7ia2 


* 

£ 


h 


a 
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These  equations  for  the  P  and  <Xare  used  to  find  the  force  required  for  the  depth 


of  the  indenter  at  the  centerline  and  the  stress  along  the  bottom  edge  of  the  indenter, 
respectively. 

2.2.2  Nonlinear  Spherical  Elastic  Indentation  Theory 

As  mentioned  in  Chapter  1,  Lin  et  al.  [12]  developed  a  finite  element  model  in 
order  to  derive  force-displacement  relations  for  spherical  indentation  of  soft  modulus 
materials  for  several  hyperelastic  functions.  For  the  present  study,  the  single-term  Ogden 
function  Lin  used  in  his  derivation  is  the  primary  focus.  The  single-term  Ogden  has  an 
energy  function  (W)  of  form: 


(2.12)  W  =  -9-(4“  +  r  +  Tz“  -  3) 

a 

(2.13)  cr=^(r-l+/ i~a'2-1) 

a 

where  po  and  a  are  fitting  parameters  and  X=  Xy=  "k7=  Xx  are  the  stretch  ratios,  go  also  has 
the  physical  meaning  of  the  initial  shear  modulus.  The  stretch  ratio  is  related  to  the  strain 
by  the  equation:  X  =  1  +  s  .  Taking  these  functions,  Lin  implemented  the  following 
approach: 

1 .  Assumed  stress  of  form 


(2.14) 


<r  =  f{C,X) 


2.  Resolved  the  sign  differences  between  standard  engineering  and  common 
indentation  notation  by  redefining 


(2.15) 


X  =  \-e 

cr*=-f(B,  l-e) 
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3.  Divide  (2. 1 5)  by  e  *  to  obtain 


(2.16) 


(2.17) 


CT  f(B,  1-g) 

£  *  £•* 

4.  For  incompressible  materials  Go  is  equal  to  Eo/3  and  taking  e*  to  0  and 
comparing  that  to  (2.8)  gives 

20  En 


3«-(l -v2) 


=  BjB 


(2.18) 


5.  Applying  the  (2.7)  yields 

P 


=  1-0.2-) 
Tra’  i? 


(2.19) 


6.  Finally,  the  contact  radius  no  longer  maintains  the  relationship  a=\l~Rh  for 
hyperelastic  functions,  therefore  a  new  relation  is  assumed 

a  =  Rx~hyhz 

where  x,  y,  and  z  are  constants.  This  relationship  was  formed  by  performing  finite 
element  analysis  of  the  scenario.  These  models  were  studied  for  insight  to  be  used  in 
developing  the  models  in  this  present  study. 

From  these  steps,  a  new  force-displacement  function  was  developed  as  described 
in  Lin,  et  a/.’s  paper; 


(2.20) 


(2.21) 


P  =  — - [(1-0.2  - )-“/2-‘  -  (1  -  0.2 -)“-*  ] 
a  R  R 


s=40£« 
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This  force-displacement  equation  is  compared  against  the  output  from  the  nonlinear 
elastic  finite  element  models. 
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2.2.3  Standard  Sharp  Indentation  Theory 


In  order  to  fully  understand  the  indentation  process,  the  theory  and  application  to 
the  elastic -plastic  materials  was  investigated.  The  finite  element  model  is  a  power-law 
hardening,  isotropic  material  that  behaves  according  to  the  following  equation: 


Where:  s  is  the  total  strain,  cy  is  the  yield  stress  of  the  material,  and  n  is  a  work  hardening 
exponent. 


Figure  15:  Elastic-plastic  Power-  Hardening  Stress  Strain  Curve 
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In  the  present  research,  a  total  of  2200  points  are  used  to  describe  the  relationship 
between  the  flow  stresses  and  the  plastic  strains  with  the  plastic  strains  within  the  range 
of  0  <  8P  <  200%.  According  to  Y.  Cao  et  al.,  [30]  this  amount  of  the  points  is  sufficient 
to  well  determine  the  plastic  behavior  of  power  law  materials. 

The  standard  method  of  using  instrumented  indentation  was  developed  by  Oliver 
and  Pharr  [3].  The  load  displacement  curve  shown  in  Chapter  1  is  reprinted  here  for 
convenience  in  Figure  16. 


> 

p  ! 


Figure  16:  Standard  Indentation  Curve 

One  of  the  most  common  indenter  geometries  is  the  Berkovich  Pyramid,  shown  in  Figure 
17.  It  is  a  three-sided  pyramid  with.  Many  indentation  experiments  are  carried  out  with  a 
Berkovich  indenter  made  of  diamond,  which  is  the  most  common  material  for  an  indenter 
tip  due  to  its  high  modulus  (1170  GPa)  and  low  Poisson  ratio  (0.07).  In  order  to  simplify 
the  analysis,  the  Berkovich  indenter  head  can  be  modeled  as  an  analytically  rigid  cone 
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with  a  half  apex  angle  of  70.3  degrees  so  that  the  cross  sectional  contact  area  is  the  same 
for  a  given  depth. 


indent  shape 


side  view 


Figure  17:  Geometry  of  a  common  Indenter  head:  Berkovich  Pyramid 

The  data  from  the  load-displacement  curve  is  used  to  calculate  the  material  properties  of 
the  specimen.  This  section  begins  with  the  equations  for  the  desired  properties  (hardness 
H  and  elastic  modulus  E)  and  then  presents  the  supporting  calculations  for  those 
equations.  Hardness  is  found  by: 

(2.23)  H=- 

A 

where  P  is  the  applied  force  and  A  is  the  projected  contact  area  (defined  in  Eqn.  11).  The 
Young’s  modulus  (the  modulus  of  elasticity,  E)  of  the  specimen  is  calculated  from  the 
reduced  modulus  Er  by  rearranging  equation  (2.2)  to  the  following; 

(2.24)  Es=(  \-vs2) 

The  reduced  modulus  can  be  found  using  the  relationship  developed  by  Oliver  and  Pharr: 
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(2.25) 


2  pjA 


Er  = 


where  S  is  the  initial  slope  of  the  unloading  curve,  also  referred  to  as  the  contact  stiffness, 
and  p  is  the  slope  of  the  indenter  tip. 

2.2.4  Elastic-Plastic  Spherical  Theory 

Furthering  the  investigation  into  standard  engineering  materials,  another 
algorithm  was  studied  and  reproduced.  This  finite  element  model  and  accompanying 
coding  analysis  is  based  on  a  system  developed  by  the  civil  engineering  department  at 
Columbia  University.  [5]  For  other  similar  processes  using  a  sharp  indenter,  in  order  to 
find  the  complete  range  of  the  stress  strain  curve  (to  include  the  plastic  regime)  several 
indentations  must  be  made  with  various  indenter  angles.  This  is  because  different 
materials  can  produce  the  same  indentation  curve.  These  are  known  as  meta-materials. 
This  can  be  very  cumbersome  and  time  consuming.  This  study  used  the  mathematical 
method  of  representative  stress  and  strain  to  find  a  unique  solution  for  the  stress-strain 
curve  with  a  single,  deep  indentation.  The  authors,  Zhao,  et  al.  defined  the  representative 

strain  to  be  the  plastic  strain,  £  . 


(2.26) 


£  =  £e  +  £  =  £e  +  £r  ( £e  is  the  elastic  strain) 


Correspondingly,  the  representative  stress  is: 


(2.27) 


CJ ,  -tL  n 


The  representative  stress  is  shown  in  Figure  15  in  the  previous  as  or  and  there  are  two 
selections  of  or  in  the  analysis.  For  spherical  indentation,  dimensional  analysis  leads  to 
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(2.28) 


C 


and 

(2.29) 


=  /(■ 


R,<?r(£r) 


,n) 


S  h  E„ 

7 TT  =  S(~, 7- -T.w) 

hER  R  <jR  \£r i 


From  equations  (2.28)  and(2.29),  Zhao  e/  al.  were  able  to  produce  the  constants  required 
to  model  the  stress-strain  curve  for  a  material  up  to  the  ultimate  yield  point.  They  used 
extensive  finite  element  analysis  to  find  equations  for  the  f  and  g.  Two  points  were 
chosen  from  the  loading  portion  of  the  force-indentation  curve:  one  at  hi  =  0. 13*h/R  and 
one  at  h2  =  0.3*h/R.  The  point  I12  also  corresponded  to  the  maximum  indentation  depth. 
Figure  3  shows  a  rough  approximation  of  their  locations.  These  2  points  were  substituted 
into  equation  (2.28)  to  produce,  along  with  (2.29),  the  3  surfaces  required  to  find  the 
solution  for  Er,  n,  and  ay. 


(2.30) 

(2.31) 


C.  =  % 

^r(^r)  2(Jr{£r) 

C  P 

(Jr(sr)  h22aR2(sR2) 


=  -/i( ~ TT  i\-w) 

°^R  \SR  / 


=  /2C 


ed 


2  V  _  2  /  2 


('.2> 


-«) 


(2.32) 


S 

h2ER 


,n) 


The  forward  analysis  used  by  Zhao,  et  al.,  produced  fittings  for  the  values  of  fl,  f2,  and 
g.  They  are: 

(2.33)  ^(m1,«)  =  /t1(^)x^(n) 
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(2.34) 


f2(m2,n)=h2(m2)xk2(n) 


(2.35) 


g(m2,n)  =  A2(n)  + 


A  i  (n)  —  A  2  (ft) 

1 +[m2  /  q(n)]p{n) 


Where  m{  = 


and  m2  = -  R 

2  _  2 


(2.36) 


hx{mx)  =  32.77  -  52.59  ln(iWj)  +  33.46(ln(w,))2 
-4.800(ln(m1  ))3  +  0.2 147(ln(iw1  ))4 


(2.37) 

(2.38) 

(2.39) 

(2.40) 

(2.41) 

(2.42) 

(2.43) 


h2(m2)  =  3.817  -12.731n(m2)  +  11.99(ln(m2))2 
-2.032(ln(m2))3  +  0.1049(ln(m2))4 

k,  {n)  =  1 .00 1  ■ -  0.26 1  On -  0.52 1  In1  +  0. 1 547«3 

k2  {n)  =  1 .002  -  0.7637/?  - 1 .920 n2  + 1 .255  ft3 

4(ft) =3.66556+0.0244179ft 

4(ft)= 6.06122-2. 15891« 

?(n)  =  29.0856- 24. 3547« 
p(n)  =  1.31861-0. 154675k 


For  this  present  study  the  reverse  analysis  Zhao,  et  al.  developed  using  these  set  of 
equations,  (2.30)  -  (2.43),  is  investigated.  This  is  discussed  further  in  Chapter  5. 


2.3  Muscle  Structure  and  Tensile  Tests  Theory 

In  order  to  understand  the  strength  of  a  muscle,  it  is  important  to  understand  how 
it  works.  The  force  required  to  stretch  a  muscle  to  a  given  length  can  be  divided  into  two 
components:  active  and  passive.  The  passive  component  is  the  contribution  from  the 
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material  structure  and  the  properties  of  the  myofilaments.  The  active  force  component  is 
the  contribution  from  the  contraction  of  the  myofilaments.  An  example  force-length 
curve  is  shown  in  Figure  18.  The  active  component  cannot  be  measured  directly.  The 
passive  and  total  forces  of  the  muscle  are  measured  separately  and  the  difference  between 
the  two  represents  the  contribution  of  the  active  state.  For  the  purposes  of  this  study,  the 
properties  of  the  individual  muscle  motor  units  are  examined  in  the  passive  state  (i.e.  no 
electrically  stimulated  contraction  of  the  muscles). 
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Figure  18:  Force-Length  Example  Diagram  for  Muscles  showing  the  Total  force- 
length  response  is  a  sum  of  the  Active  and  Passive  Properties  of  the  muscle 

Passive  muscle  tissue  is  a  viscoelastic,  hyperelastic,  anisotropic  material. 
However,  after  preconditioning,  the  viscoelastic  nature  of  the  material  becomes  minimal 
and  the  material  is  then  regarded  as  pseudo-elastic.  Pseudo-elastic  materials  behave  as 
one  elastic  material  during  loading  and  another  one  during  unloading.  [5]  Additionally, 
the  organization  of  the  muscle  fibers  simplifies  the  material  to  a  transversely  isotropic 
material  and  is  assumed  to  be  incompressible. 
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In  addition  to  the  indentation  tests,  the  theory  of  which  has  already  been  discussed 
at  length,  uniaxial  tensile  tests  were  also  attempted  on  the  muscle  samples.  Tensile  tests 
involve  applying  a  force  to  a  sample  in  one  direction  and  observing  the  change  in  its 
length.  It  is  one  of  the  fundamental  tests  of  materials  strength.  For  a  uniaxial  load  the 
material  described  in  equation  (2.1)  simplifies  the  engineering  stress  and  engineering 
strain  described  by: 

(2.44)  a  =  Ee 

A L  P 

where  £  =  —  and  cr= — .  Lo  and  Ao  are  the  initial  length  and  area  of  the  specimen 

4  A 

being  tested  [31]. 

2.4  Summary 

In  this  Chapter  the  mathematical  formulations  for  the  finite  element  model 
materials  were  outlined.  The  derivations  of  the  nonlinear  elastic  and  the  elastic-plastic 
models  were  discussed.  These  equations  and  the  Hertzian  linear  elastic  equations  are 
used  to  verify  the  models  developed  in  Chapter  3. 

In  addition  to  the  mathematical  formulas  for  indentation,  the  structure  of  the 
hawkmoth  muscle  was  explained.  The  flight  mechanics  of  the  moth  and  how  the  muscles 
drive  that  motion  were  diagramed.  This  information  is  used  to  analyze  the  experimental 
data. 
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III.  Finite  Element  Model  Development 


3.1  Chapter  Overview 

In  this  chapter  the  development  and  validation  of  the  FE  models  are  discussed. 
Issues  such  as  element  type  and  number,  boundary  conditions,  analysis  techniques,  and 
contact  type  are  analyzed.  The  mesh  resolution  of  the  sample  and  the  given  indenter  for 
the  elastic  sample  will  be  compared  to  Hertz  contact  equations  from  Chapter  2  to  show 
convergence. 

3.2  Finite  Element  Overview 

In  this  section  the  finite  element  models  of  the  three  indenter  head  geometries 
probing  into  an  elastic  half-space  is  explained.  The  factors  going  into  the  mesh  are 
explored  and  the  results  are  validated  against  the  analytic  Hertzian  solutions  to  ensure 
mesh  refinement  so  that  the  models  may  be  used  in  the  analysis  section  in  Chapter  5. 

The  commercial  finite  element  analysis  software  package  ABAQUS,  version  6.10 
is  used  in  the  indentation  simulation.  Information  in  this  section  about  the  finite  element 
model  comes  from  the  ABAQUS  Users  Manual  [32],  ABAQUS  is  a  suite  of  powerful 
engineering  simulation  programs,  based  on  the  finite  element  method,  which  can  solve 
problems  ranging  from  relatively  simple  linear  analysis  to  the  most  challenging  nonlinear 
simulations. 

3.3  Analysis  Considerations 

ABAQUS  consists  of  two  main  analysis  modules:  ABAQUS/Standard  and 

ABAQUS/Explicit.  ABAQUS/Standard  is  a  general-purpose  analysis  module  that  can 

solve  a  wide  range  of  linear  and  nonlinear  problems  efficiently,  accurately  and  reliably. 
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ABAQUS/Explicit  is  a  special-purpose  analysis  module  that  uses  an  explicit  dynamic 
finite  element  formulation.  It  is  suitable  for  short,  transient  dynamic  events.  The 
indentation  procedure  is  assumed  to  be  quasi-static  problem,  in  which  no  rate  effect  is 
considered,  therefore  ABAQUS/Standard  is  employed  in  this  work. 

In  the  indentation  simulation,  there  are  two  sources  of  nonlinearity:  material 
nonlinearity  and  geometric  nonlinearity.  The  indentation  procedure  can  produce  large 
deformation  in  the  solids  underneath  and  near  the  indenter.  The  magnitude  of 
displacement  can  affect  the  response  of  the  structure  (geometry  nonlinearity).  ABAQUS 
uses  the  modified  Newton-Raphson  method  to  obtain  solutions  for  nonlinear  problems  in 
the  Standard  module. 

3.4  Models 

For  all  of  the  simulations  in  this  study  the  indenter  head  is  much  stiffer  than  the 
medium  being  indented.  This  allows  the  indenter  head  to  be  modeled  as  an  analytically 
rigid  solid.  The  rigid  surface  is  associated  with  a  rigid  body  reference  node,  whose 
motion  governs  the  motion  of  the  surface.  Since  only  one  node  is  computed,  this  saves 
computer  resources  and  simplifies  analysis. 

Both  the  indenter  and  half  space  are  modeled  as  axisymmetric  geometry.  The  bottom 
of  the  model  is  fixed  and  the  remaining  two  sides  are  free  as  shown  in  Figure  19.  Four-node 
axisymmetric  linear  quadrilateral  elements  are  utilized  for  the  half-space.  Reduced 
integration  is  employed  to  spare  calculation  time.  The  element  type  for  the  elastic  and 
elastic-plastic  materials  used  in  ABAQUS  is  ‘CAX4R’,  in  which  the  letter  or  number 
indicates  the  element  is  continuum,  axisymmetric,  4-node  bilinear,  reduced  integration 
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with  hourglass  control  respectively.  This  element  is  not  able  to  be  used  for  hyperelastic 
material  problems.  The  hybrid  version  of  the  CAX4R  is  used  in  this  case  called 
CAX4RH.  The  axisymmetric  elements  are  shown  in  Figure  20. 


Figure  19:  Model  Boundary  Conditions:  Y-axis  Symmetry  (Left  side)  and  Fixed 

(Bottom) 
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4  3 


1  2 


Figure  20:  Diagram  of  simplification  made  by  assuming  no  variation  in  the  angular 
direction  for  the  CAX4R  and  CAX4RH  elements 

The  mesh  is  designed  so  that  the  meshing  is  refined  near  the  indenter  (in  order  to 
resolve  the  contact  conditions  and  allow  for  accurate  contact  area  determination),  but  also 
is  sufficiently  large  so  that  it  approximates  a  semi-infinite  solid.  Accordingly,  the  mesh  is 
chosen  large  enough  for  each  calculation  so  that  the  results  obtained  are  insensitive  to  the 
movement  of  the  outer  boundaries  of  the  mesh.  For  all  three  models,  a  structured  grid  is 
used  in  order  to  decrease  computational  time  and  to  resolve  the  fine  mesh  required  at  the 
point  of  indentation.  The  grid  can  be  seen  in  Figure  21  and  Figure  22  for  the  flat  punch 
and  sphere,  respectively.  Size  of  the  sample  half  space  (LxW)  in  the  figures)  is 
determined  for  each  scenario  due  to  varying  indentation  depths  and  material  effects. 
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Figure  21:  ABAQUS  Axisymmetric  Model  of  Cylindrical  Flat  Punch  indenting  a 
sample  with  partition  lines  to  divide  sample  in  order  to  refine  the  mesh.  Inset: 

axisymmetric  assumption  diagram 
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Figure  22:  ABAQUS  Axisymmetric  Model  of  Sphere  indenting  a  sample  with 
partition  lines  to  divide  sample  in  order  to  refine  the  mesh.  Inset:  axisymmetric 

assumption  diagram 


A  3D  representation  of  the  2  types  of  models  can  be  seen  in  Figure  23.  An  example  mesh 
can  be  seen  in  Figure  24  with  note  of  the  additional  refinement  of  the  mesh  on  the  right- 
hand  side  of  the  figure. 


Figure  23:  Close  up  View  of  3D  revolution  of  Axisymmetric  Models: 
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Flat  Punch  (Left)  &  Spherical  (Right) 


Figure  24:  Example  mesh  of  a  sample  as  a  whole  (left)  and  close  up  showing  the 
refined  area  underneath  the  indenter  (right) 


3.7  Summary 

In  summary,  for  a  given  scenario,  an  axisymmetric  model  is  used  with  either  a 
spherical  or  cylindrical  flat  punch  indenter.  The  scenarios  are  considered  to  be  quasi¬ 
static  to  simulate  the  rate-independent  nature  of  preconditioned  biological  tissues.  The 
mesh  and  dimensions  of  the  models  are  varied  according  to  the  experiment  to  be 
evaluated. 
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IV.  Experimentation  Methodology 


4.1  Chapter  Overview 

In  this  Chapter  the  methodology  for  the  attempted  experiments  is  outlined.  The 
thought  behind  this  methodology  is  explained  and  the  assumptions  for  the  data  gathered 
are  discussed. 

4.2  Nanoindentation  Experimentation 

Hawkmoth  pupae  were  maintained  in  a  temperature  controlled  room  with  a  14 
hour  on/  10  hour  off  light/dark  cycle  until  they  hatched.  More  information  of  the  raising 
of  the  hawkmoth  can  be  found  in  the  appendix.  Once  the  moths  hatched,  they  were 
moved  to  a  secondary  cage  until  they  were  needed.  The  moths  could  last  up  to 
approximately  a  week  in  this  cage  before  they  died. 

Prior  to  the  dissection,  a  glass  slide  was  attached  to  a  metal  puck  with 
Crystalbond™  (Crystalbond,  AREMCO,  Valley  Cottage,  NY).  (Figure  25) 
Crystalbond™  is  a  heat-activated  adhesive  material.  A  thin  ring  was  attached  to  the  glass 
slide  with  cyanoacrylate  adhesive  (super  glue).  This  ring  contains  the  saline  solution  so 
that  the  sample  will  not  desiccate  during  the  experiment.  Small  pieces  of  glass  are  also 
adhered  to  glass  slide.  These  are  approximately  the  same  height  as  the  tissue  sample.  It 
was  hoped  that  these  glass  blocks  could  be  used  as  a  guide  for  the  nanoindenter  to  gauge 
the  distance  to  the  surface. 
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Figure  25:  Sample  Puck  for  indentation  experiment  with  rubber  ring  attached  to 

glass  slide  to  hold  saline  solution 

The  moths  were  taken  out  of  their  cage  approximately  1  hour  before  the  test  was 
conducted  and  asphyxiated  using  a  paper  towel  dipped  in  acetone  in  a  small  closed  off 
container.  This  was  to  maintain  the  freshness  of  the  sample  and  an  attempt  to  prevent 
decomposition  as  much  as  possible. 

Once  the  insects  were  dead,  the  wings,  head,  abdomen,  and  legs  were  all  removed 
using  small  scissors.  This  left  only  the  thorax  remaining.  The  thorax  was  held  under  a 
slow  stream  of  running  water  and  scrubbed  with  a  toothbrush  to  remove  the  small  scales. 
This  allowed  the  exoskeleton  to  be  clearly  visible  to  aid  in  the  dissection.  The  scissors 
were  again  used  to  clip  the  thorax  along  the  dashed  line  shown  in  Figure  26. 
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Figure  26:  Partially  Dissected  Moth  highlighting  the  point  of  incision  of  the 

exoskeleton  (wings  are  removed) 


The  bottom  of  the  exoskeleton  was  removed,  leaving  the  top  of  the  exoskeleton 
and  exposed  muscles.  A  cut  was  made  at  the  front  and  rear  of  the  thorax,  slicing  through 
the  DLMs.  As  the  DLMs  are  not  attached  to  the  top  of  the  thorax,  this  allowed  the 
bundle  to  be  removed  with  tweezers.  Care  was  taken  to  cut  the  muscles  as  close  to  the 
exoskeleton  as  possible  to  remove  the  largest  amount  of  tissue.  Isolated  muscle  tissue  can 
be  seen  in  Figure  27.  Drops  of  saline  solution  were  applied  to  the  tissue  samples  to 
maintain  moisture  until  they  were  ready. 
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Figure  27:  Isolated  Individual  muscle  units  of  the  Hawkmoth  DLMs 

When  ready  to  test,  the  muscle  sample  was  patted  dry  to  remove  excess  moisture 
and  weighed  to  the  nearest  0.001  g.  A  small  bit  of  cyanoacrylate  was  placed  on  the  glass 
slide  next  to  glass  surface  find  aide.  The  muscle  sample  was  placed  on  top  of  the 
adhesive  as  shown  in  Figure  28.  The  sample  was  covered  with  a  moist  paper  towel  and 
allowed  to  sit  for  several  minutes.  This  allowed  the  adhesive  to  dry  while  maintaining 
the  moisture  of  the  muscle.  Once  dry,  the  ring  was  filled  with  saline.  This  setup  is  much 
like  that  used  by  V.T.  Nayar  et  al.  [13]  indentation  of  porcine  sclera  mentioned  in 
Chapter  1 . 
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Figure  28:  (left)  Cartoon  showing  sample/adhesive  interaction  and  (right)  muscle 
sample  in  solution  adhered  to  slide  with  glass  block  flush  with  its  right  edge. 


The  puck  as  a  whole  was  weighed  in  order  to  be  weighed  again  following  the  test 
to  determine  if  there  was  any  water  lost  from  the  sample.  The  sample  puck  was  then 
placed  in  the  Nanoindenter  G200  test  rack  as  shown  in  Figure  22. 


Figure  29:  (a)  Agilent  G200  Nanoindenter  (b)  Sample  Puck  in  Test  Rack  of  Indenter 
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Two  indenter  heads  were  used:  120pm  and  300pm  radius  spherical  head.  The 
Standard  XP  Indentation  Head  was  used.  The  surface  sensitivity  was  set  at  “very  fine” 
with  an  approach  distance  of  1  mm.  Through  the  nanovideo  camera  microscope, 
locations  on  the  sample  and  on  the  glass  near  the  sample  were  chosen.  The  indenter  was 
programmed  to  indent  the  glass  and  then  use  the  height  location  to  zero  in  on  the  surface 
of  the  sample.  It  was  unable  to  find  the  surface  of  the  muscle  for  either  the  120pm  or 
300pm  indenter  head.  Discussion  of  this  outcome  is  in  Chapter  5. 

4.3  Uniaxial  Tensile  Experimentation 

The  second  experiment  with  the  hawkmoth  muscle  fiber  examined  the  tensile 
properties  of  the  specimen.  Dissection  and  isolation  of  the  muscle  units  occurred  as  in 
section  3.3.1.  Each  sample  was  weighed  and  measured.  Length  and  the  width  at  three 
discrete  points  were  measured  using  a  digital  microscope.  The  three  widths  were 
averaged  and  treated  as  approximately  the  diameter  muscle.  This  diameter  was  used  to 
find  the  cross  sectional  area  used  for  the  stress  calculations. 

The  material  was  tested  using  the  Agilent  T150  Universal  Testing  Machine 
(UTM)  shown  in  Figure  30  along  with  its  specifications.  The  T150  has  the  necessary 
displacement  and  load  resolution  required  for  a  sample  of  this  size  with  this  low  of  an 
expected  modulus. 
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Agilent  T150  UTM  Specifications 

Maximum  load 

500mN  (50.8gm) 

Load  resolution 

50nN  (5.1  pgm) 

Maximum  actuating  transducer  displacement 

±1mm 

Displacement  resolution 

<0.1  nm 

Dynamic  displacement  resolution 

<0.001nm 

Maximum  crosshead  extension 

200  mm 

Extension  resolution 

35nm 

Extension  rate 

0.5pm/s  to  5mm/s 

Dynamic  frequency  range  (sample  dependent) 

0.1  Hz  to  2.5kHz 

CDA  Option 

Force  amplitude  range 

0.1  pN  to  4.5mN 

Frequency  range  characterization  of  instrument 
dynamic  response  (sample  dependent) 

0.01Hz  to  200  Hz 

Figure  30:  Agilent  T150  UTM  and  Specifications 

In  order  to  attach  the  samples  to  the  machine,  an  interface  is  needed  between  the 
grips  and  the  tissue  sample.  As  was  done  by  Calvo,  et  al.  [17],  the  specimen  is  attached 
using  cyanoacrylate  to  small  pieces  of  sandpaper  which  in  turn  are  attached  to  paper 
templates.  These  pieces  of  sandpaper  prevent  slippage  between  the  sample  and  the 
clamps.  The  applicability  of  the  sandpaper/glue  combination  is  confirmed  in  Ng  et  al. 
study  in  2005  [33].  The  paper  templates  allow  for  the  handling  of  the  samples  with 
placing  too  much  stress  on  them  before  the  testing  starts.  The  templates  are  connected 
with  tape.  Once  the  sample  is  clamped  into  the  testing  apparatus,  the  tape  is  cut,  allowing 
for  the  machine  stress-strain  analysis  to  start.  Setup  can  be  seen  in  Figure  31. 
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Figure  31:  Template  card  technique  showing  Sample  in  Testing  Machine.  Edges  of 
card  have  been  cut,  allowing  the  sample  to  be  stretched. 


Due  to  the  size  limitations  of  the  sample  and  the  restrictive  nature  of  the  template, 
each  of  the  gauge  lengths  for  the  tests  was  restricted  to  approximately  2  mm.  Sample 
hydration  was  limited  to  the  beginning  of  each,  much  in  the  way  that  Gras,  et  al.  [16],  did 
in  their  testing  of  the  stemocleidomastoideus  muscle.  This  lack  of  hydration  throughout 
the  test  resulted  in  the  inability  to  perform  a  preconditioning  cycle  to  the  specimen. 
According  to  theory  in  Chapter  2,  the  elastic  properties  of  the  muscle  would  vary  with 
different  strain  rates  due  to  the  viscoelastic  nature  of  the  unconditioned  material.  To 
investigate  this,  three  strain  rates  were  to  be  tested  (3e-4,  le-4,  and  le-3  mm/mm/s), 
however,  the  uniaxial  testing  was  only  completed  for  4  samples  at  strain  rates  of  3e-4 
mm/mm/s.  The  clamps  for  the  testing  machine  split  and  new  ones  were  unable  to  be 
procured  in  time  to  complete  the  remainder  of  the  planned  tests.  Results  that  were 
obtained  as  well  as  improvements  on  the  testing  method  are  discussed  in  Chapter  5. 
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V.  Analysis  and  Results 


5.1  Chapter  Overview 

In  this  Chapter,  the  validity  and  usefulness  of  the  finite  element  models  are 
analyzed  and  their  results  commented  on.  The  nanoindentation  experiment  was  analyzed 
to  determine  why  the  experiment  failed  and  whether  the  Agilent  G200  is  useful  for 
materials  with  similar  properties  as  the  muscle.  Results  from  the  uniaxial  tension  tests 
are  reported  and  the  results  critiqued  to  determine  better  testing  methods. 


5.3  Results  of  Finite  Element  Simulations 

In  this  section,  3  scenarios  are  analyzed.  The  first  looks  at  the  algorithm 
developed  for  elastic -plastic  spherical  indentations.  The  second  looks  at  the  porcine 
sclera  experiment  discussed  in  the  literature  review  and  how  boundary  effects  could  come 
into  play  for  similar  future  experiments.  The  third  example  looks  at  the  effect  of  friction 
on  the  results  of  Lin,  et  al.  in  their  analysis  of  mouse  cartilage. 

5.3.1  Spherical  Indenter  into  Elastic  Plastic  Medium 

In  this  scenario,  the  elastic  plastic  properties  of  a  material  are  determined  by  the 
analysis  method  developed  by  Zhao  et  al.  [5]  described  in  Chapter  2.  The  forward 
analysis  fitting  functions  shown  in  equations  (2.30)  -  (2.43)  are  used.  Equations  (2.30) 
through  (2.32)  are  repeated  here  for  convenience. 
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The  reverse  analysis  requires  inputs  for  Pi,  P2,  hi,  I12,  and  S.  These  inputs  can  be 
measured  from  an  experiment.  However,  in  the  absence  of  experimental  data  finite 
element  analysis  may  be  used  to  produce  these  values.  Materials  with  known  values  for 
Gy,  E,  and  n  (u  is  assumed  to  be  0.3)  can  be  input  into  a  finite  element  model.  The  force- 
indentation  data  may  be  extracted  from  the  results  of  this  analysis  and  the  values  for  Pi, 
P2,  hi,  h2,  and  S  can  be  found. 

For  example,  the  properties  for  A533-B  steel  are  approximately  E  =  210  GPa,  ay  = 
400  MPa,  and  n  =  0.127  [5],  These  values  were  input  into  the  model  outlined  in  Chapter 
3  for  the  spherical  indentation.  The  radius  of  the  indenter  was  set  to  be  788  pm  and  the 
sample  to  be  50  times  the  radius  to  remove  any  boundary  condition  effects. 

The  first  step  was  to  check  the  mesh  refinement  against  the  spherical  elastic 
analytic  equations  shown  in  Chapter  2.  The  depth  of  indentation  is  set  to  the  maximum 
value  required  for  the  algorithm,  h=0.3R,  or  h  =  236.4  pm.  This  indentation  can  be  seen 
in  Figure  32. 
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Figure  32:  Spherical  Elastic  Maximum  Indentation,  2D  (left)  and  3D  (right).  Colors 
indicate  stress  in  the  vertical  direction  (MPa) 


The  force  indentation  curve  for  4  mesh  sizes  are  compared  in  Figure  33.  The  ‘100x100’ 
refers  to  the  number  of  elements  in  the  comer  square  of  the  mesh.  The  larger  the  number 
is,  the  finer  the  mesh. 


Figure  33:  Mesh  Resolution  Comparison,  F-d  Elastic  Spherical 
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The  stress  under  the  indenter  at  the  maximum  depth  is  compared  with  the  analytic 
solution  in  Figure  34.  The  maximum  stress  at  the  middle  of  the  indenter  is  predicted 
better  using  the  coarse  mesh.  Due  to  there  being  such  a  large  displacement,  the  material  is 
out  of  the  linear  range.  Since  the  force-displacement  diagram  is  the  driver  for  the 
analysis,  that  refinement  is  the  primary  mesh. 


Figure  34:  Mesh  Resolution  Comparison,  F-d  Elastic  Spherical 

The  finest  mesh  was  considered  to  have  converged;  therefore  it  was  used  to  determine  the 
values  required  for  the  analysis.  The  material  properties  of  the  model  were  changed  to 
include  the  work  hardening  strain  values  (E=210  GPa,  ay,=400MPa,  and  n  =  0.127)  and 
the  indentation  analysis  was  run  again.  The  stress  field  can  be  seen  in  at  the  maximum 
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indentation  depth  in  Figure  35  and  after  the  indenter  has  separated  from  the  surface  in 
Figure  36.  Residual  stresses  can  be  seen  in  the  in  the  sample.  The  force  indentation  curve 
for  this  scenario  can  be  seen  in  Figure  37. 


Figure  35:  Spherical  Plastic  Maximum  Indentation,  2D  (left)  and  3D  (right) 


Figure  36:  Spherical  Plastic  Residual  Indentation,  2D  (left)  and  3D  (right) 
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Figure  37:  Force  Displacement  Curve  for  Forward  Analysis  (R=788  pm) 

The  values  for  Pi,  P2,  hi,  h2,  and  S  are  extracted  from  Figure  37  and  are  shown  in  Table 
1 .  These  values  are  used  to  determine  the  constants  to  be  input  into  the  algorithm  shown 
in  Table  2. 


Table  1:  Measured  Values  from  Force  Indentation  Curve 


hi  (pm) 

h2  (pm) 

PI  (N) 

P2(N) 

102.44 

236.4 

1023.8 

2280.4 

Table  2:  Input  Values  to  Zhao  Algorithm 


Cl 

(GPa) 

C2 

(GPa) 

S 

(MPa*m) 

97.56 

40.81 

349.3 
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From  here  the  flow  chart  shown  in  Figure  38  may  be  followed  to  determine  the 
material  properties.  The  steps  in  this  flow  chart  are  solved  using  the  code  found  in 
Appendix  B.  The  range  for  E/cy  was  initially  set  to  be  100  to  1000,  the  range  for  n  was  0 
to  0.6,  and  the  ay/ao  range  was  set  to  0.5  to  1.5.  This  allowed  for  a  wide  range  of 
engineering  materials  to  be  in  consideration  for  this  analysis. 


Figure  38:  Flow  Chart  for  Determining  Material  Properties 

The  output  values  from  the  analysis  are  listed  in  Figure  39.  In  Table  3,  the  results 
for  each  of  the  materials  are  compared  to  those  found  by  Zhao,  et  al.  in  their  model. 
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Figure  39:  Comparison  of  Force  Indentation  for  Forward  and  Reverse  Analysis 
Table  3:  Zhao  Algorithm  Output  Comparison  with  Symbolic  Solver 


E  (GPa) 

oy  (MPa) 

n 

Input  Parameter 

210 

400 

0.127 

Zhao  Output 

206 

402 

0.125 

Percent 

Difference 

2% 

1% 

2% 

Dauby  Output 

247 

396 

0.125 

Percent 

Difference 

18% 

1% 

2% 

Although,  he  was  only  off  by  2%  for  this  material,  Zhao  reported  values  with 
consistently  10%  error  of  for  each  value.  The  first  analysis  of  the  present  study  had 
slightly  higher  errors  for  E,  consistently  8  or  so  %  higher  than  Zhao.  One  possible  reason 
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for  the  difference  could  be  the  choice  of  a  symbolic  solver.  During  the  analysis,  in  order 
to  find  the  reference  stress  a  nonlinear,  explicit  equation  must  be  solved.  The  analysis 
was  then  repeated  using  a  numerical  solver  to  determine  the  reference  stress.  These 
results  are  shown  in  Table  4.  The  E  value  produced  was  more  in  line  with  input  values 
with  slight  uptick  in  the  yield  stress  and  n  value. 

Table  4:  Zhao  Algorithm  Output  Comparison  with  Numerical  Solver 


E  (GPa) 

oy  (MPa) 

n 

Input  Parameter 

210 

400 

0.127 

Zhao  Output 

206 

402 

0.125 

Percent 

Difference 

2% 

1% 

2% 

Dauby  Output 

235 

380 

0.13 

Percent 

Difference 

12% 

2% 

5% 

One  potential  application  for  this  technique  to  biological  tissues  is  for  harder 
materials  such  bone.  The  hawkmoth  exoskeleton  would  seem  to  be  a  possible  candidate 
for  this  analysis.  Since  the  approximate  thickness  of  the  exoskeleton  in  30  pm  [21],  there 
would  be  limitations  to  the  size  of  the  indenter  when  setting  up  this  experiment.  In  order 
to  model  a  half  space  approximately  50  times  the  maximum  depth  of  the  indentation 
(0.3R),  one  would  be  limited  to  an  indenter  with  a  radius  of  approximately  2  pm.  Any 
larger  and  substrate  effects  would  be  at  risk  of  influencing  results. 

5.3.2  Flat  Punch  to  Elastic  Medium 

In  this  scenario,  how  the  boundary  affects  the  experiment  conducted  by  Nayar  et 
al.  [13]  was  analyzed  using  the  elastic  finite  element  model.  As  mentioned  previously, 
this  experiment  involved  an  80  micron  diameter  cylindrical  flat  punch  indenting  porcine 
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sclera.  The  indentation  chosen  was  only  4  pm  deep,  therefore,  according  to  (2.1 1),  the 
strain  resulting  from  this  experiment  is  approximately  0. 1  and  falls  within  the  accepted 
region  for  elasticity.  The  linearity  of  this  region  could  be  verified  by  a  hyperelastic  finite 
element  analysis;  however,  the  only  data  that  could  be  found  for  the  sclera  was  elastic. 

In  their  experiment,  the  sclera  samples  were  cut  into  1  cm  squares  with  a 
thickness  of  approximately  1  to  1 .2  mm.  This  would  appear  to  be  a  large  enough  sample 
to  remove  effects  of  the  boundary  on  the  analysis.  However,  finite  element  model  can 
help  in  ensuring  that  there  is  no  discrepancy.  For  the  purpose  of  this  simulation,  the 
sample  was  set  to  1.2  mm  (1200  pm)  on  each  side  to  match  the  thickness  of  the  real  life 
sample. 

First,  the  flat  punch  model  was  modified  to  simulate  the  scenario  of  the  sclera 
experiment  so  that  the  geometry  (40  pm  radius  indenter)  and  material  properties 
(E=30kPa)  are  made  to  match.  The  mesh  is  varied  and  the  force  displacement  curve  and 
the  stress  field  under  the  indenter  are  mapped  in  Figure  40  and  Figure  41,  respectively. 
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Figure  40:  Mesh  Refinement  for  Flat  Punch  using  Force-Displacement  Relation 


Figure  41:  Mesh  Refinement  for  Flat  Punch  using  Stress-Radius  Relation 
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The  force-displacement  curve  produced  stiffer  results  as  the  mesh  became  finer, 
however  the  stress  field  under  the  indenter  more  closely  matched  the  analytic  solution. 
The  stress  field  is  the  primary  means  of  measuring  the  boundary  effects  so  the  finer  mesh 
was  used  as  the  model. 

Next  the  boundary  conditions  were  modified  to  have  rollers  on  the  left  and  right 
side  of  the  sample  as  in  Figure  42.  This  modification  verifies  that  the  1200  by  1200  pm 
sample  adequately  represents  an  infinite  half-space.  Additionally,  the  sample  was 
increased  to  4800  by  4800  pm  as  a  second  check.  These  results  are  used  for  comparison 
as  the  width  and  depth  of  the  sample  are  adjusted. 
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Figure  42:  Model  Boundary  Conditions:  Y-axis  Symmetry  (Left  side,  Right  side) 

and  Fixed  (bottom) 

As  shown  in  Figure  21,  Figure  22,  and  Figure  42,  W  represents  the  horizontal 
distance  from  the  point  of  indention  to  the  vertical  boundary  and  L  represents  the  vertical 
distance  from  the  surface  of  the  sample  to  the  bottom  of  the  sample.  The  bottom  is  fixed, 
representing  a  hard  surface  underneath  the  sample.  Typically  the  rule  of  thumb  for 
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thickness  of  a  sample  is  to  maintain  an  indentation  depth  that  is  no  more  than  10-20%  of 
the  thickness  of  the  sample.  [34] 

To  analyze  what  effect  the  horizontal  boundary  distance  has  on  the  measurement 
of  the  elastic  modulus,  the  simulation  was  run  using  W  values  of  600  pm  and  300  pm  and 
again  with  L  values  of  600  pm  and  300  pm.  The  stress  underneath  the  indenter  and  the 
force-displacement  relationship  were  plotted  in  Figure  43  and  Figure  44,  respectively,  for 
varying  W  values.  The  stress  and  force-displacement  relationship  for  varying  L  values 
are  shown  in  Figure  45 and  Figure  46. 


Figure  43:  Comparison  of  the  stress  along  the  bottom  of  a  flat  punch  indenting 

several  sample  sizes 
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Figure  44:  Comparison  of  the  force-displacement  relationship  for  a  flat  punch 
indenting  several  sample  sizes  with  varying  horizontal  distance  (W) 


Figure  45:  Comparison  of  the  stress  along  the  bottom  of  a  flat  punch  indenting 

several  sample  sizes 
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Figure  46:  Comparison  of  the  force-displacement  relationship  for  a  flat  punch 
indenting  several  sample  sizes  with  varying  vertical  distance  (L) 

In  Table  5,  the  effects  of  the  smallest  boundary  condition  are  compared  to  the 
large  sample  and  the  Hertzian  solution.  The  Hertzian  solution  is  treated  as  truth  for  the 
percent  error  calculations.  As  can  be  seen  from  Table  5,  the  smaller  boundary  can  have  a 
large  effect  on  the  calculation  of  the  elastic  modulus.  Errors  were  most  noticeable  for  the 
thickness  of  the  sample  at  21.9%.  This  shows  that  for  low  modulus  materials,  the 
thickness  of  the  sample  should  be  more  carefully  monitored  than  for  a  similarly  sized 
stiffer  material. 


65 


Table  5:  Summary  of  effects  of  E  from  the  smallest  sample  boundary  length 
compared  to  Hertzian  analytic  solution  and  a  large  sample 


a (kPa) 
(Midpoint) 

%  Error 

E (kPa) 

%  Error 

Hertz 

-1.3929 

0.0% 

30.0 

0.0% 

1200x1200 

-1.3144 

5.6% 

30.1 

0.4% 

1200x300 

-1.5065 

8.2% 

34.2 

14.1% 

300x1200 

-1.0368 

25.6% 

23.4 

21.9% 

5.3.3  Spherical  Indenter  into  Hyperelastic  Medium 

In  this  scenario,  a  rigid  spherical  indenter  is  indented  into  a  hyperelastic  half- 
space.  Up  until  now,  the  coefficient  of  friction  between  the  sample  and  the  indenter  has 
been  assumed  to  be  frictionless.  For  the  vast  majority  of  indentation  experiments  the 
effect  of  friction  is  assumed  to  be  negligible.  According  to  Shacham  S,  et  al.,  [35]  for 
muscle  tissues  in  contact  with  bone  the  coefficient  of  friction  could  be  as  high  as  0.36. 
This  is  not  a  one-to-one  relationship  with  a  diamond  or  sapphire  indenter  in  contact  with 
a  tissue,  however,  the  relationship  this  factor  has  for  a  hyperelastic  material  is  worth 
investigating. 

The  model  begins  with  the  frictionless  scenario  established  previously.  The  mesh 
from  section  5.3.1  that  was  verified  using  the  elastic  space  is  used  for  this  analysis  as 
well.  The  indenter  size  is  adjusted  to  5  pm  and  the  sample  size  was  correspondingly 
decreased  to  125  pm  square.  Per  Lin  et  al.  experiment,  the  indentation  is  1  pm. 

In  order  to  gauge  the  effects  of  friction  on  this  model,  the  coefficient  was 
increased  to  0.1,  0.2  and  0.3.  The  indentation  stress  field  is  shown  in  Figure  47. 
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Figure  47:  Stress  Field  for  Hyperelastic  Material 

The  effects  on  the  force  indentation  curve  are  shown  in  Figure  48.  The  elastic  force- 
displacement  model  and  the  hertz  analytic  solution  were  also  plotted  for  comparison. 
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Figure  48:  Friction  Comparison  Force  Indentation  Curve  (h=lpm) 

As  can  be  seen  from  Figure  48,  friction  had  no  effect  on  the  force-displacement 
relation  whatsoever.  In  Figure  49,  the  contact  area  is  compared  for  the  differing  friction 
values  as  well.  Again,  there  is  no  change.  The  frictionless  assumption  is  valid  for  this 
indentation  depth. 
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Figure  49:  Friction  Comparison  Contact  Radius  (h=lpm) 

In  order  to  see  how  the  material  behaves  outside  of  the  indentation  depth  made  by  Lin,  et 
al,  the  scenario  was  conducted  again  for  an  indentation  of  4  pm.  As  can  be  seen  in 
Figure  50  and  Figure  51,  there  is  some  minor  variation  in  the  force  indentation  and  the 
contact  radius  from  the  finite  element  simulation.  It  appears  that  the  frictionless 
assumption  becomes  marginally  less  valid  as  the  indentation  depth  increases. 
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Figure  50:  Friction  Comparison  Force  Indentation  (h=4pm) 


Figure  51:  Friction  Comparison  Contact  Radius  (h=4pm) 
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5.3.4  Extraction  of  Ogden  Parameters  from  Force-Indentation  Curve 

The  spherical  indentation  model  described  n  section  5.3.3  can  also  be  used  to 
extract  parameters  from  experimental  data.  As  in  section  5.3.1,  experimental  data  was 
not  available  so  the  finite  element  model  was  used  to  produce  simulated  data.  Initial 
shear  modulus  (po)  and  curvature  parameters  (a)  were  input  in  the  ABAQUS  model.  The 
parameters  chose  were:  po  =  7.89kPa  and  a  =  20.  The  radius  of  the  spherical  indenter 
was  set  at  R  =  120  pm  and  size  of  the  sample  was  set  at  50  times  the  radius.  The 
indentation  depth  was  set  at  0.4R,  or  48  pm.  The  resulting  stress-strain  curve  can  be  seen 
in  Figure  52  and  long  with  the  solution  for  the  linear  elastic  medium  shown  in  blue  for 
comparison.  The  nonlinearity  of  the  solution  is  apparent. 


Figure  52:  Representative  stress  and  strain  from  a  spherical  indenter  into  a 
hyperelastic  medium  simulation  compared  against  an  elastic  medium 
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In  order  to  determine  the  coefficients  or  this  simulated  data,  the  force- 


displacement  relationship  in  equation  (2.20)  presented  by  Lin,  et  al.  was  fitted  to  the  data 


using  a  nonlinear-least  squares  data  fit.  The  code  for  this  extraction  is  in  Appendix  B. 


(2.20)  P  =  ^-[(1  -  0.2  — )"“/2“1  _  (l  _  0.2  —  )“_1  ] 

a  R  R 


Equation  (2.20)  was  rearranged  to  express  the  stress-strain  relationship  using  the 


equations  a*  =  P  /  na 2  and  s *  =  0.2  *  a  /  R  to  the  form: 


(5.1)  o-*  =  -  [(1  -e  )-“/2_1  -  (1  -  e  )a~x  ] 

a 


The  fit  to  this  equation  is  shown  in  Figure  53. 


Figure  53:  Representative  stress  and  strain  from  a  spherical  indenter  into  a 
hyperelastic  medium  simulation  with  a  curve  fit  to  extract  Ogden  parameters 


72 


The  parameters  as  determined  by  the  least  squares  method  are:  po  =  7.369  kPa  and  a  = 
23.2.  This  represents  a  7.55%  and  16.10%  percent  difference  from  the  input  po  and  a, 
respectively.  The  new  parameters  were  input  into  the  finite  element  model  to  compare 
against  the  initial  simulation.  The  results  are  shown  in  Figure  54  and  match  well. 


Figure  54;  Comparison  of  the  simulation  results  from  initial  material  parameters 
and  material  parameters  from  a  curve  fit 

5.4  Results  of  Nanoindenter  Experiment 

No  data  were  able  to  be  gathered  from  the  nanoindentation  experimentation.  For 
each  indentation  test,  the  indenter  head  would  approach  the  surface  and  push  through. 
However,  the  indenter  instrumentation  would  be  unable  to  record  that  the  surface  had 
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made  contact.  It  soon  became  apparent  the  primary  reason  was  that  the  load  vs. 
displacement  slope  for  the  material  was  too  low.  The  surface  stiffness  of  the  sample  was 
too  low  to  escape  the  noise  of  the  machine. 


Load  vs  Disp  Slope  (NAn) 


Figure  55:  Noise  from  Load  vs  Displacement  Channel  for  G200  during  Surface  Find 

To  investigate  whether  this  was  the  cause,  the  elastic  Hertz  contact  equations 
were  revisited.  Since  the  surface  contact  involves  the  initial  depth  of  indentation,  the 
material  can  be  assumed  to  be  elastic  without  the  hyperelastic  effects  dominating. 
Considering  first  the  force-displacement  equation  (reprinted  from  Chapter  2): 

(2.5)  p  =  —ErRP1P 

And  taking  the  derivative  with  respect  to  depth  yields 

dP  11 

(5.2)  lI=2E«Rh 

oh 
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Substituting  in  the  definition  of  reduced  modulus,  ER  =£7(1— v2)  andu=0.5 
(incompressible  assumption)  gives; 

dP  8Z J.’ 

(5.3)  — =-ER2h 2 

V  7  dh  3 

This  function  was  plotted  by  varying  the  elastic  modulus  while  holding  the  indenter 
radius  constant  for  select  values.  This  can  be  seen  in  Figure  56.  Additionally,  from 
Figure  55,  the  approximate  noise  level  for  the  Agilent  G200  was  approximately  10  N/m 
and  also  was  plotted  as  the  thick  horizontal  line. 


Figure  56:  Load/Displacement  Slope  vs.  Elastic  Modulus  for  Varied  Indenter 
Radius  with  noise  level  shown  as  black  horizontal  line. 
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From  the  graph,  the  smallest  modulus  that  would  be  able  to  be  measured  for  this  range  of 
indenter  sizes  would  be  approximately  3  to  6  MPa.  This  is  approximately  100  times  the 
values  for  soft  tissue  found  in  the  literature. 

Returning  to  the  hyperelastic  spherical  indentation  model  from  section  5.3.4,  the 
force  indentation  relationship  can  be  seen  in  Figure  57.The  initial  slope  for  the  first 
0.5pm  of  the  indentation  is  approximately  0.3.  This  shows  that  the  soft  material  used  in 
the  Lin,  et  al.  experiment  would  not  be  able  to  be  measured  with  the  G200  nanoindenter. 
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Figure  57:  Force-displacement  relationship  for  spherical  indenting  a  hyperelastic 
medium  with  1st  order  Ogden  potential  parameters  pO=7.97kPa  and  a=20 
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Improvements  to  the  equipment  for  this  test  are  needed.  One  possibility  is  to  use 
the  DCM II  Indentation  Head  Option  for  the  Agilent  G200.  This  option  provides  greater 
resolution  for  the  indentations  and  sensitivity  in  the  analysis  with  load  resolution  of  3  nN. 
However,  this  option  was  not  functioning  for  the  test  system  and  future  researchers  would 
need  to  have  it  replaced. 

The  more  likely  solution  to  this  problem  is  to  use  AFM  instrumentation  to 
conduct  the  experiment.  AFM  instruments  have  a  much  lower  force  and  displacement 
resolution.  Spherical  tips  are  available  for  AFM  machines  as  well,  so  the  finite  element 
development  for  the  indentation  experiments  would  still  be  applicable,  only  on  a  smaller 
scale. 

5.5  Results  of  Tensile  Experiment 

The  results  of  the  tensile  experimental  investigation  are  presented  in  this  section. 
As  mentioned  at  the  end  of  Chapter  4,  the  clamps  for  the  T 150  broke  before  the 
completion  of  the  planned  tests  could  be  completed.  A  summary  of  those  tests  that  were 
able  to  be  completed  is  given  in  Table  6.  It  should  be  noted  that  Test  1  was  aborted  after 
a  strain  of  only  approximately  0.03  due  to  a  disturbance  to  the  test  fixture.  The  results 
have  also  been  presented  in  Figure  58.  The  elastic  modulus  was  extracted  by  applying  a 
linear  bet  fit  line  to  the  portion  of  data.  The  overall  modulus  was  calculated  by  applying 
linear  regression  to  the  entire  data  set  (only  the  linear  portions  of  the  data  were  used). 

Table  6:  Results  for  4  Tests  of  Initial  Modulus  Elastic  for  Hawkmoth  Muscle 


Test 

E,  kPa 

A0,mmA2 

10,  mm 

Test  1 

802.17 

1.5394 

1.6 

Test  2 

39.80 

1.131 

1.6 
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Test  3 

998.32 

1.3273 

2 

Test  4 

801.92 

1.5394 

2 

Figure  58:  Stress-Strain  Diagram  for  all  4  tests  of  Hawkmoth  DLM  motor  unit  and 
their  individual  linear  regression  fits  and  the  overall  regression  fit. 

As  can  be  seen  from  the  Table  6  and  from  Figure  58,  3  of  the  tests  (1,3,  and  4) 
had  reasonably  agreeable  values  of  elastic  modulus  in  the  high  100s  of  kPa.  The  other 
test  (2)  is  approximately  20  times  smaller  than  the  other  three  tests.  Looking  at  Table  7, 
overall  average  of  the  slopes  and  the  standard  deviation  (SD)  of  the  slopes  is  reported. 
The  overall  regression  uses  all  the  points  in  the  sample  If  Test  2  (T2)  is  removed  because 
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it  falls  outside  of  standard  deviation  of  the  remaining  tests,  the  overall  elastic  modulus 
becomes  982  kPa.  These  results  are  plotted  in  Figure  59. 


Table  7:  Summary  Statistics  of  the  4  Uniaxial  Tension  Tests 


Test 

E,  kPa 

A0,mmA2 

10,  mm 

Overall  Average 

660.56 

1.38 

1.80 

Overall  SD 

424.05 

0.20 

0.23 

Overall  Regression 

343.72 

Average  w/o  T2 

867.47 

1.47 

1.87 

SD  without  T2 

113.32 

0.12 

0.23 

Regression  w/o  T2 

982.3 

Figure  59:  Hawkmoth  DLM  Stress-strain  curve  with  Outlier  Removed 
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However,  the  Test  2  outlier  may  not  necessarily  have  been  an  invalid  test.  Most 
importantly,  the  small  sample  size  precludes  making  too  many  inferences  about  which 
test  is  valid  and  which  is  not.  Additionally,  the  level  of  hydration  of  the  samples  can 
cause  variability  in  their  elastic  properties.  A  desiccation  test  was  devised  to  determine 
how  much  weight  the  muscles  lost  when  exposed  to  the  air.  A  single  muscle  unit  was 
extracted  as  described  in  Section  4.2.  It  was  patted  dry  to  remove  excess  moisture  and 
placed  on  a  scale  and  weighed  to  the  nearest  milligram.  The  sample  was  reweighed  at 
regular  intervals  for  a  period  of  16  minutes.  The  results  of  this  experiment  can  be  seen  in 
Figure  60.  The  muscle  unit  lost  over  half  its  weight  in  only  a  period  of  16  minutes.  On 
average  the  tests  lasted  approximately  6  minutes  there  for  approximately  12-15%  of  the 
muscles  water  weight  could  have  evaporated  during  that  time. 


80 


8 


U) 

E 


O) 

'q3 


£ 


0  2  4  6  8  10  12  14  16 

time,  min 


Figure  60:  Chart  of  desiccation  of  Hawkmoth  muscle  over  time  while  exposed  to  air 

Qualitatively,  as  the  samples  become  less  hydrated,  they  become  much  harder  to 
the  touch  and  would  most  likely  have  a  higher  value  for  E.  It  is  possible  that  the  other  3 
tests  lost  more  moisture  than  test  2.  Additionally,  this  idea  that  the  other  three  lost 
moisture  looks  like  more  of  a  possibility  when  compared  to  the  results  of  other  tensile 
tests  of  muscles  as  shown  in  Table  8. 
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Table  8:  Summary  Comparison  with  prior  Muscle  Tensile  Experiments 


Study 

Modulus 

(kPa) 

Animal 

Gras ,etal.  [16] 

111 

Human 

Dorfman  et  al.  [18] 

2.34 

Caterpillar 

Collinsworth,  et  al. 

[37] 

45.6 

Rabbit 

Fung  [6] 

20-160 

Rabbit 

Dauby 

334 

Moth 

Better  ways  of  keeping  the  sample  hydrated  need  to  be  used  for  this  type  of  test. 
The  sample  dried  out  prior  to  the  conclusion  of  the  test.  As  mentioned  previously,  the 
hydration  method  used  was  similar  to  that  used  by  Gras,  et  al.  Their  study  was  able  to 
hydrate  only  at  certain  points  in  the  test  because  their  samples  were  much  larger.  The 
smaller  surface  area  to  volume  ratio  of  their  samples  did  not  allow  for  as  much 
desiccation.  A  better  way  of  testing  the  samples  would  be  to  keep  them  completely 
submerged  in  the  saline  solution. 

The  small  length  of  the  sample  with  respect  to  the  diameter  can  also  affect  the 
results.  In  order  to  estimate  what  this  effect  might  be,  a  new  finite  element  code  was 
developed.  This  model  was  again  axisymmetric,  with  element  type  CA4XR.  The  first 
1.5  mm  of  the  sample  was  fixed  along  the  edges  of  one  end  while  the  displacement  was 
applied  to  the  edges  of  the  other  end.  This  simulates  the  glue  holding  the  ends  of  the 
sample.  The  mesh  was  refined  near  these  boundary  points.  The  mesh  and  boundary 
conditions  for  this  model  can  be  seen  in  Figure  61.  The  diameter  of  each  sample  was 
fixed  at  2mm,  which  is  approximately  the  diameter  of  the  muscle  unit.  Four  different 
gauge  lengths  were  analyzed:  2,  4,  7,  and  10  mm. 
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Figure  61:  Boundary  Conditions  (left)  and  Mesh  (right)  for  Uniaxial  Tension  Model 

(2mm  gauge  length  shown) 

In  the  simulations  in  Figure  62,  the  bottom  right  figure  has  a  gauge  length  of  2mm 
which  is  approximately  the  situation  in  this  experiment.  For  a  large  strain  of  0.5,  there  is 
a  much  larger  stress  variation  across  the  midsection  of  the  sample.  The  gauge  length  is 
increased  going  clockwise.  As  can  be  seen,  the  further  away  from  the  boundary 
conditions,  the  more  uniform  the  stress  field  in  the  sample. 
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Figure  62:  Stress  cross  sections  for  all  tests  to  50%  strain.  (Clockwise  from  bottom 
left:  2mm  gauge  length,  4  mm  gauge  length,  7  mm  gauge  length,  and  10  mm  gauge 

length.) 

New  ways  of  securing  the  sample  are  recommended  by  this  report  to  be 
investigated.  One  possible  way  is  to  dissect  the  moth  in  a  way  to  keep  the  attachment 
points  of  the  exoskeleton  intact.  This  would  keep  the  full  10mm  length  of  the  muscle  to 
be  used  as  the  gauge  length.  Also  the  boundary  conditions  would  remain  intact,  thereby 
better  simulating  the  natural  environment.  Better  dissecting  skills  would  be  required  to 
accomplish  this  new  method. 


5.6  Summary 

Three  finite  element  models  were  developed  with  possible  applications  to 
biological  tissue  and  possible  applications  to  the  hawkmoth.  One  model  determined  the 
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elastic-plastic  properties  of  a  wide  range  of  materials  from  one  single  spherical 
indentation.  The  second  examined  the  boundary  effects  of  the  experiment  with  porcine 
sclera.  The  third  examined  the  effect  of  coefficient  of  friction  on  a  hyperelastic  material. 

The  experimentation  portion  of  the  analysis  returned  a  mixed  bag  of  results.  The 
nanoindentation  experiment  was  unable  to  gather  any  data,  although  an  upper  bound  was 
found  on  the  indentation  modulus  of  a  material  able  to  be  characterized  by  the  Agilent 
G200  Nanoindenter  with  a  spherical  tip.  The  tensile  testing  was  able  to  gather  data  on 
the  longitudinal  modulus  of  the  hawkmoth  muscle.  Compared  to  literature,  the  modulus 
was  high  most  likely  due  to  desiccation  of  the  sample  on  the  test  device. 

VI.  Conclusions  and  Recommendations 


6.1  Chapter  Overview 

In  this  chapter  the  conclusions  that  were  obtained  as  a  result  of  this  research  are 
outlined.  Additionally,  recommendations  for  future  research  are  discussed.  A  summary 
of  the  study  concludes  the  report. 

6.2  Conclusions  of  Research 

A  finite  element  model  was  developed  to  analyze  the  elastic,  power  law¬ 
hardening  properties  of  a  wide  range  of  material  properties.  The  values  reported  from  the 
analysis  differed  from  the  actual  values  by  approximately  18%.  This  was  slightly  higher 
than  the  reported  values  from  Zhao,  et  al. 

The  second  finite  element  model  was  developed  to  analyze  the  boundary  effects 
of  an  experiment  analyzing  the  flat  punch  indentation  of  the  material.  It  showed  that  the 
sample  Nayer,  et  al.  used  was  more  than  satisfactory  for  conducting  the  experiment. 
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Additionally  it  was  shown  that  reducing  the  horizontal  boundary  to  300  pm  could  cause 
as  much  as  a  23%  error  in  the  evaluation  of  the  elastic  modulus.  The  third  model  looked 
in  the  experiment  by  Lin,  et  al.  to  see  what  effects  friction  may  have  had  on  the  analysis. 
For  the  indentation  depth  in  the  Lin  study,  friction  would  have  had  little  to  no  impact  on 
the  results.  Had  the  choice  of  indentation  been  deeper,  a  more  pronounced  effect  could 
have  been  seen. 

The  experimental  nanoindentation  experiment  was  unable  to  gather  data  due  to 
the  limitations  in  the  instrument  in  measuring  a  material  with  as  low  a  modulus  of 
elasticity  as  the  hawkmoth.  An  upper  bound  on  the  modulus  was  established  of 
approximately  3  MPa.  This  is  consistent  with  literature  values  of  modulus  in  the  range  of 
1-50  kPa  for  soft  tissues  from  other  experimentation. 

The  uniaxial  tension  test  was  able  to  map  the  stress  strain  curve  for  strains  up  to 
10%.  Reported  initial  modulus  of  elasticity  values  were  343  kPa  although  only  three 
experiments  and  one  partial  experiment  were  able  to  be  conducted  before  the  test  fixture 
broke.  Additional  concerns  for  the  testing  methodology  are  centered  on  the  dehydration 
of  the  samples  during  testing.  This  loss  of  water  most  likely  had  an  increase  in  the 
stiffness  of  the  samples  and  larger  vales  for  modulus.  An  additional  concern  for  the  test 
involves  the  small  gauge  length  with  respect  to  size  of  the  specimen.  Finite  element 
models  show  probable  uneven  stress  values  throughout  the  specimen. 

6.3  Recommendations  for  Future  Research 

The  algorithm  developed  for  the  elastic-plastic  model  could  be  applied  to  a 
variety  of  materials.  However,  in  relation  to  the  problem  of  the  material  properties  of  the 
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hawkmoth,  it  could  be  used  to  determine  the  entire  stress  strain  curve  if  instrumentation 
could  be  found  with  an  indenter  of  radius  small  enough  to  go  0.3R  into  the  material  while 
not  going  into  the  substrate  of  the  moth  exoskeleton. 

For  future  research  involving  indentation  experiments  with  biological  tissues,  the 
Agilent  G200  with  the  Standard  XP  Indentation  Head  should  not  be  used.  It  should  not  be 
used  because  the  load-displacement  slope  of  the  material  is  too  low  to  be  recorded  by  the 
machine.  Other  options  such  as  an  Atomic  Force  Microscope  should  be  considered. 
Should  the  Nanoindentation  experimentation  data  become  available,  the  finite  element 
models  should  be  used  to  analyze  and  isolate  the  material  properties. 

Future  research  into  the  uniaxial  tension  test  should  account  the  hydration  of  the 
sample  better.  Applying  saline  solution  directly  prior  to  the  test  is  not  adequate  and 
spraying  the  samples  during  test  causes  the  test  to  be  aborted.  A  horizontal  testing 
apparatus  with  the  sample  completely  submerged  would  be  a  preferable  solution  to  the 
current  testing  system.  This  improved  hydration  would  allow  for  preconditioning  of  the 
sample.  Additionally,  improved  dissection  technique  could  allow  for  a  longer  gauge 
length  which  would  improve  both  handling  of  the  samples  and  reduce  boundary  condition 
effects. 

6.4  Summary 

In  summary,  this  study  described  a  software  package  utilizing  the  commercial 
finite  element  suite  ABAQUS  to  allow  hyperelastic  materials  to  be  considered  with  the 
application  towards  soft  biological  tissue. 
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Two  soft  tissue  models  were  developed  and  used  to  investigate  boundary  effects 
and  coefficient  of  friction  on  two  experiments  conducted  by  outside  researchers.  The 
third  model  used  a  technique  by  Zhao,  et  al.  to  determine  the  elastic-plastic  properties  of 
Another  model  successfully  computed  the  material  properties  of  an  elastic -plastic 
material  using  only  one  spherical  indentation.  This  model  could  be  applied  to  the 
hawkmoth  exoskeleton. 

Two  experiments  were  attempted.  The  first,  a  Nanoindentation  experiment  with 
the  flight  muscle  of  a  hawkmoth  was  unsuccessful.  The  instrumentation  was  unable  to 
measure  the  modulus  of  the  material.  The  second  experiment  was  a  uniaxial  tension  test 
on  the  muscle.  This  experiment  was  able  to  obtain  an  initial  elastic  modus  of  the 
material;  however  the  results  may  skew  high  due  to  loss  of  moisture  in  the  sample  during 
the  experiment. 
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Appendix  A:  Finite  Element  ABAQUS  Input  Files 
Sphere  (Elastic-Plastic) 

^Heading 

**  Job  name:  zhao788  Model  name:  zhao-788mm 
**  Generated  by:  Abaqus/CAE  6.10-1 

^Preprint,  echo=NO,  model=NO,  history=NO,  contact=NO 
*  * 

**  PARTS 
*  * 

*Part ,  name=Sample 
*End  Part 

■k  * 

*Part ,  name=Spherical 

*End  Part 
*  * 

*  * 

**  ASSEMBLY 
*  * 

"'"Assembly,  name=Assembly 

~k  ~k 

* Instance,  name=Sample-l ,  part=Sample 

0.,  -40246.7628172297,  0. 

*Node 


1, 

1970 . , 

38276.7617 

2, 

1970 . , 

40246.7617 

3, 

o., 

40246.7617 

4, 

o., 

38276.7617 

5, 

39400. , 

38276.7617 

6, 

39400 . , 

40246.7617 

etc 

**  INTERACTION  PROPERTIES 

■k  k 

"^Surface  Interaction,  name=IntProp-l 

1., 

*Friction,  slip  tolerance=0 . 005 

0., 

* Surface  Behavior,  pres  sure  -over  closure=HARD 

k  k 

**  BOUNDARY  CONDITIONS 

k  k 

**  Name:  Axi  Type:  Symmetry/Antisymmetry/Encastre 
"^Boundary 

_PickedSet33 ,  YASYMM 

**  Name:  Bottom  Type:  Symmetry/Antisymmetry/Encastre 
"^Boundary 

_PickedSet34 ,  ENCASTRE 
*  * 

**  INTERACTIONS 
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*  * 


**  Interaction:  Int-1 

^Contact  Pair,  interaction=IntProp-l ,  type=SURFACE  TO 
SURFACE 

Sample-1 . "Sample  top".  Spherical . "Sphere  Surf" 

*  *  _ 


*  * 

**  STEP:  Down 
*  * 

*Step,  name=Down,  nlgeom=YES,  inc=10000 
^Static 

0.001,  1.,  le-06,  0.01 
*  * 

**  BOUNDARY  CONDITIONS 
*  * 

**  Name:  Down  Type:  Displacement/Rotation 
^Boundary 

_PickedSet21,  1,  1 
_PickedSet21,  2,  2,  -236.4 

_PickedSet21,  6,  6 
*  * 

**  OUTPUT  REQUESTS 
*  * 

^Restart,  write,  frequency=0 
*  * 

**  FIELD  OUTPUT:  F-Output-1 
*  * 

^Output,  field,  variable=PRESELECT 
*  * 

**  HISTORY  OUTPUT:  Down 
*  * 

^Output,  history 

*Node  Output,  nset=Spherical . "Sphere  RP  Set" 

RF2,  U2 

*Node  Print,  nset=Spherical . "Sphere  RP  Set",  SUMMARY=NO 
RF2,  U2 
*End  Step 

*  *  _ 


*  * 

**  STEP:  Up 
*  * 

*Step,  name=Up,  nlgeom=YES,  inc=10000 
^Static 

0.001,  1.,  le-05 ,  0.05 
*  * 

**  BOUNDARY  CONDITIONS 
*  * 

**  Name:  Axi  Type:  Symmetry/Antisymmetry/Encastre 
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^Boundary,  op=NEW 
_PickedSet33 ,  YASYMM 

**  Name:  Bottom  Type:  Symmetry/Antisymmetry/Encastre 
^Boundary,  op=NEW 
_PickedSet34 ,  ENCASTRE 

**  Name:  Down  Type:  Displacement/Rotation 
*Boundary,  op=NEW 

**  Name:  Up  Type:  Displacement/Rotation 
*Boundary,  op=NEW 
_PickedSet22 ,  1,  1 
_PickedSet22,  2,  2,  236.4 

_PickedSet22 ,  6,  6 

*  * 

**  OUTPUT  REQUESTS 
*  * 

*Restart,  write,  frequency=0 
*  * 

**  FIELD  OUTPUT:  F-Output-1 
*  * 

*Output ,  field,  variable=PRE SELECT 
*  * 

**  HISTORY  OUTPUT:  Down,  Up 
*  * 

^Output,  history 

*Node  Output,  nset=Spherical . "Sphere  RP  Set" 

RF2,  U2 

*Node  Print,  nset=Spherical . "Sphere  RP  Set",  SUMMARY=NO 
RF2,  U2 
*End  Step 

Flat  Punch  (Elastic) 

^Heading 

**  Job  name:  FP4-ref inelOO  Model  name:  Flat  Punch-ref ine- 
Copy 

**  Generated  by:  Abaqus/CAE  6.10-2 

^Preprint,  echo=NO,  model=NO,  history=NO,  contact=NO 
*  * 

**  PARTS 
*  * 

*Part,  name=Cone 

*End  Part 
*  * 

*Part,  name=Sample 

*End  Part 
*  * 

*  * 

**  ASSEMBLY 
*  * 

^Assembly,  name=Assembly 
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*  * 


^Instance,  name=Cone-l,  part=Cone 

0., 

-1500 . 

r 

0. 

*Node 

1, 

o., 

1500  .  , 

0. 

*Nset,  nset=Cone-l-Ref Pt  , 

internal 

i/ 

*Nset,  nset="Cone  RP  Set" 

1, 

^Surface,  type 

=SEGMENTS ,  name=MCone 

surf" 

START, 

40., 

1580  . 

LINE, 

40., 

1501 . 

CIRCL, 

39., 

1500  .  , 

39., 

1501  . 

LINE, 

o., 

1500  . 

*Rigid  Body,  ref  node=Cone 

:-l-Ref Pt  , 

analytical 

surf ace="Cone 

surf" 

^End  Instance 

*  * 

* Instance,  name=Sample-l , 

part=Sample 

0., 

-3846.76281722965, 

0. 

*Node 

1, 

o 

00 

3766.7627 

2, 

o 

00 

3846.7627 

3, 

0., 

3846.7627 

4, 

0., 

3766.7627 

5, 

800  .  , 

3766.7627 

6, 

800.  , 

3846.7627 

v. 

o 

00 

3046.7627 

8, 

0., 

3046.7627 

9, 

800.  , 

3046.7627 

More  nodes,  etc 

203,  204,  205, 

206,  207 

*Elset,  elset= 

surface  set. 

instance= 

:Sample-l ,  generate 

100,  10000 

,  100 

*End  Assembly 
*  * 

**  MATERIALS 
*  * 

*Material,  name=Sample 
^Density 
2 . 65e-15 , 7 0  . 

^Elastic 

0.03,  0.499 
*  * 

**  INTERACTION  PROPERTIES 
*  * 

* Surface  Interaction,  name=IntProp-l 

1 .  , 
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*Friction,  slip  tolerance=0 . 005 

0.1, 

^Surface  Behavior,  pressure-overclosure=HARD 
*  * 

**  BOUNDARY  CONDITIONS 
*  * 

**  Name:  Axi  Type:  Symmetry/Antisymmetry/Encastre 
^Boundary 

_PickedSet8 ,  YASYMM 

**  Name:  Bottom  Type:  Symmetry/Antisymmetry/Encastre 
^Boundary 

_PickedSet7 ,  ENCASTRE 
*  * 

**  INTERACTIONS 
*  * 

**  Interaction:  Int-1 

^Contact  Pair,  interaction=IntProp-l ,  type=SURFACE  TO 
SURFACE 

Sample-1 . "Sample  top",  Cone-1. "Cone  surf" 

*  *  _ 


*  * 

**  STEP:  Down 
*  * 

*Step,  name=Down,  nlgeom=YES,  inc=1000 
^Static 

0.01,  1.,  le-05,  0.1 
*  * 


**  BOUNDARY  CONDITIONS 
*  * 


**  Name:  Down  Type: 
^Boundary 

_PickedSet39,  1,  1 
_PickedSet39,  2,  2, 

_PickedSet39,  6,  6 
*  * 

**  OUTPUT  REQUESTS 
*  * 

^Restart,  write, 

*  * 

**  FIELD  OUTPUT: 

*  * 

^Output,  field, 

*  * 


Displacement /Rotation 

-4  . 


f requency=0 
F-Output-1 
variable=PRE SELECT 


**  HISTORY  OUTPUT:  Contactl 
*  * 

^Output,  history 
^Contact  Output 
CARE A, 

*  * 
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**  HISTORY  OUTPUT:  Down 
*  * 

*Node  Output,  nset=Cone-l . "Cone  RP  Set” 

RF2,  U2 

*Node  Print,  nset=Cone-l . "Cone  RP  Set",  summary=no 
RF2,  U2 

^Contact  Print,  summary=no 
CARE A, 

*End  Step 

*  *  _ 


*  * 

**  STEP:  Up 
*  * 

*Step,  name=Up,  nlgeom=YES,  inc=1000 
^Static 

0.01,  1.,  le-05 ,  0.1 
* 

**  BOUNDARY  CONDITIONS 
*  * 

**  Name:  Axi  Type:  Symmetry/Antisymmetry/Encastre 
^Boundary,  op=NEW 
_PickedSet8 ,  YASYMM 

**  Name:  Bottom  Type:  Symmetry/Antisymmetry/Encastre 
^Boundary,  op=NEW 
_PickedSet7 ,  ENCASTRE 

**  Name:  Down  Type:  Displacement/Rotation 
^Boundary,  op=NEW 

**  Name:  Up  Type:  Displacement/Rotation 
^Boundary,  op=NEW 
_PickedSet40,  1,  1 
_PickedSet40,  2,  2,  4. 

_PickedSet40,  6,  6 
*  * 

**  OUTPUT  REQUESTS 
*  * 

^Restart,  write,  frequency=0 
*  * 

**  FIELD  OUTPUT:  F-Output-1 
*  * 

^Output,  field,  variable=PRESELECT 
*  * 

**  HISTORY  OUTPUT:  Contactl 
*  * 

^Output,  history 
^Contact  Output 
CARE A, 

*  * 

**  HISTORY  OUTPUT:  Down 
*  * 
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*Node  Output,  nset=Cone-l . "Cone  RP  Set” 
RF2,  U2 
*End  Step 


Sphere  (Hyperelastic) 

^Heading 

**  Job  name:  Spheref inalcf 0  Model  name:  Sphere_new 
**  Generated  by:  Abaqus/CAE  6.10-2 

^Preprint,  echo=NO,  model=NO,  history=NO,  contact=NO 
*  * 

**  PARTS 
*  * 

*Part,  name=Cone 

*End  Part 
*  * 

*Part,  name=Sample 

*End  Part 
*  * 

*  * 

**  ASSEMBLY 
*  * 

^Assembly,  name=Assembly 
*  * 

^Instance,  name=Cone-l,  part=Cone 

0.,  -1500.,  0. 

*Node 

1,  -9. 18485047e-16,  1500.,  0. 

*Nset,  nset=Cone-l-Ref Pt_,  internal 

*Nset,  nset="Cone  RP  Set” 

1, 

^Surface,  type=SEGMENTS,  name="Cone  surf” 

START,  2.63692888321968,  1509.24812971375 
CIRCL ,  0.,  1500.,  0., 

1505. 

*Rigid  Body,  ref  node=Cone-l-Ref Pt_,  analytical 
surf ace=” Cone  surf" 

*End  Instance 
*  * 

* Instance,  name=Sample-l ,  part=Sample 

0.,  -3846.76281722965,  0. 

*Node 


1, 

12., 

3834 .7627 

2, 

12., 

3846.7627 

3, 

o., 

3846.7627 

4, 

o., 

3834 .7627 

5, 

125  .  , 

3834 .7627 

6, 

125  .  , 

3846.7627 

95 

7, 

12., 

3721.7627 

8, 

o., 

3721.7627 

9, 

125  .  , 

3721.7627 

10, 

12., 

3834.88281 

11, 

12., 

3835.00293 

12, 

12., 

3835.1228 

13, 

12., 

3835.24292 

14, 

12., 

3835.36279 

15, 

12., 

3835.48291 

16, 

12., 

3835.60278 

17, 

12., 

3835.7229 

18, 

12., 

3835.84277 

19, 

More  nodes,  etc 

12., 

3835.96289 

*Nset,  nset=_PickedSet3 9 ,  internal,  instance=Cone-l 

Ir 

*Nset,  nset=_PickedSet4 0 ,  internal,  instance=Cone-l 

If 

^End  Assembly 
*  * 

**  MATERIALS 
*  * 

*Material,  name=Sample 
* Hype re las tic,  ogden 

0.0143,  7.3,  0. 

*  * 

**  INTERACTION  PROPERTIES 
*  * 

^Surface  Interaction,  name=IntProp-l 

1  •  , 

^Friction,  slip  tolerance=0 . 005 

0., 

* Surface  Behavior,  pres sure -over closure=HARD 
*  * 

**  BOUNDARY  CONDITIONS 
*  * 

**  Name:  Axi  Type:  Symmetry/Antisymmetry/Encastre 
^Boundary 

_PickedSet8 ,  YASYMM 

**  Name:  Bottom  Type:  Symmetry/Antisymmetry/Encastre 
^Boundary 

_PickedSet7 ,  ENCASTRE 
*  * 

**  INTERACTIONS 
*  * 

**  Interaction:  Int-1 

^Contact  Pair,  interaction=IntProp-l ,  type=SURFACE  TO 
SURFACE 

Sample-1 . "Sample  top",  Cone-1. "Cone  surf" 
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*  * 


*  * 

**  STEP:  Down 
*  * 

*Step,  name=Down,  nlgeom=YES,  inc=1000 
^Static 

0.01,  1.,  le-05,  0.1 
*  * 

**  BOUNDARY  CONDITIONS 
*  * 

**  Name:  Down  Type:  Displacement/Rotation 
^Boundary 

_PickedSet39,  1,  1 
_PickedSet39,  2,  2,  -1. 

_PickedSet39,  6,  6 
*  * 

**  OUTPUT  REQUESTS 
*  * 

^Restart,  write,  frequency=0 
*  * 

**  FIELD  OUTPUT:  F-Output-1 
*  * 

^Output,  field,  variable=PRE SELECT 
*  * 

**  HISTORY  OUTPUT:  contact 
*  * 

^Output,  history 
^Contact  Output 
CARE A, 

*  * 

**  HISTORY  OUTPUT:  Down 
*  * 

*Node  Output,  nset=Cone-l . "Cone  RP  Set” 

RF2,  U2 

*Node  Print,  nset  =  Cone-1. "Cone  RP  Set",  summary=no 
RF2,  U2 

^Contact  Print,  summary=no 

CAREA 

*End  Step 
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Appendix  B:  MATLAB  Codes 


Spherical  (Elastic  Plastic) 

This  code  will  find  the  values  from  the  force-indentation  curve 

clear 
hold  on 

file3  =  ' zhao788elastic40 . dat ' % [ ' SphereE ' 

num2str (ctrE) ,  '  s  '  , num2str (ctrs) ,  '  n ' ,  num2str (ctrn) ,  '  . dat 1 ] ; 

f id=f open ( f ile3 ,  ' r ' ) ;  %  Opens  .dat  file  for  Reading 

ctr  =  1; 

ctr2  =  1; 

while  ~feof (fid) 

tline=fgetl (fid) ; 

if  isempty ( strf ind ( tline ,  ' N  ODE  OUTPUT'))  ==0  %Scan  till 
header 

for  i=l : 9 

tline=fgetl (fid) ;  %  Read  10  lines  of  junk  lines 

end 

dataLine=f getl ( f id) ; %Grabs  the  reaction  force  and  displacement 
data  =  sscanf (dataLine ,  ' %i  %f  %f'); 

P (ctr) =-data (2 ) /10A6; %  extracts  force  data  and  converts  to  N 
del (ctr) =-data (3) /10A6;  %  extracts  disp  data  and  converts  to  m 
ctr  =  ctr+1; 

end 

if  isempty ( strf ind (tline, 'C  ONTACT  OUTPUT'))  ==  0  %Scan 
till  header 

for  i=l:ll 

tline=f getl ( f id) ;  %  Read  10  lines  of  junk  lines 

end 

dataLine=fgetl ( f id) %Grabs  the  contact  stress 
if  dataLine  >0 

%  data  =  sscanf (dataLine,  ' %i  %s  %f  % f ' ) 
data  =  sscanf (dataLine ,  ' %f ' ) 

CAd (ctr2) =data (1) / (10A6) A2  %  extracts  contact  area  data  and 
converts  to  mA2 
end 

ctr2  =  ctr2#l; 

end 

end 

f close (fid) ; 

R  =  788; 

del  =  del (find (del>=0 ) ) ^10 A  6 ; 

P  =  P (find(del>=0) ) ; 
plot (del , P, 's') 

max_del  =  max (del); 

max_P  =  max (P) ; 

ind_max_P  =  f ind ( P==max_P) ; 

ind_zero_P  =  find(P==0); 

unload_P  =  P ( ind_max_P : ind_zero_P ( 1 ) ) ; 

unload_del  =  del ( ind_max_P : ind_zero_P ( 1 ) ) ; 
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load_P  =  P (1 : ind_max_P) ; 
load_del  =  del (1 : ind_max_P) ; 

%  figure 

%  plot (unload_del-unload_del (end) , unload_P,  'o') 

S  =  [diff (unload_P) ./diff (unload_del) ] ' 

dif f_del_13=abs ( ( 0 . 13*R-load_del ) ) ' 

ind_de 1_1 3=f ind ( di f  f _de 1_1 3==min (dif f_del_13 ) ) 

del_13=load_del (ind_del_13) 

P_13  =  load_P (ind_del_13) 

s  =  fitoptions ( 'Method' , ' NonlinearLeastSquares ' , . . . 

' Lower ',[0,0],... 

'Upper' , [max_del , max_P] , . . . 

' StartPoint ' ,  [1.1  31],... 

' TolX ' ,  eps ) 
form  =  [  ' B*dAm ' ] 

f  =  f ittype ( f orm,  'independent',  'd',  'coefficients',  {'m','B'}, 

' options ' , s )  ; 

fitobj  =  fit ( (unload_del-unload_del (end) ) ' ,  unload_P',  f) 

cf  =  coef fvalues ( fitobj ) 

max_P 

max_del 

S_real  =  le6*cf (1) *cf (2) * (max_del-unload_del (end) ) A (cf (1) -1) 
hold  on 

%  plot (unload_del-unload_del (end) , cf (2) * (unload_del- 
unload_del (end) )  . Acf  (1) ,  'k') 

%101. 1,1003  103.1,1034 

slope  *  (1034-1003) / (103.1-101.1) ; 

P_13_int  =  slope* ( 102 . 44-101 . 1 )  +1003  %  Interpolates  between  the  points 
to  get  at  0 . 13R 


This  code  is  for  solving  the  Zhao  algorithm 

%Starting  values 
%  %fine9 

E_sigy  =  linspace ( 300 , 800 , 100 ) ; 
n  =  . 1 : .005: .15 

sigy_sig0  =  linspace (. 5 , 1 . 5 , 100 ) ; 
sigO  =  400e6; 

v  =  0.3; 

eRl  =  0.0374; 
eR2  =  0.0674; 

ER_sigy  =  E_sigy/ ( l-vA2 ) ; 

R  =  7  8  8  e  -  6 ; _ 
dell  =  0 . 1 3  *  R ; 
del 2  =  0 . 3  *  R ; 

%  Outputs  from  FE  analysis  or  from  experiment 
PI  =  1 . 0238e+03  +0  ; 
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P2  =  2 . 2804e+03  +0  ; 

S  =  3 . 4 932e+08  +0  ; 

Cl  =  P 1 / de 1 1 A  2  ; 

C2  =  P2/del2A2; 


%  make  1  if  it  is  the  first  time  running  the  program 
%  make  it  any  number  greater  than  1  after  that 
run_num  =  12 

%Builds  the  reference  stress  matrix  and  saves  to  matrix. mat 
if  run_num==l 

sigRl  =  ones ( length (E_sigy) , length ( sigy_sig0 ) , length (n) ) ; 
sigR2  =  ones ( length (E_sigy) , length ( sigy_sig0 ), length (n) )  ; 
Ep  =  ones ( length (E_sigy) , length ( sigy_sig0 ), length (n) ) ; 

Sp  =  ones ( length (E_sigy) , length ( sigy_sig0 ), length (n) ) / 
np  =  ones ( length (E_sigy) , length ( sigy_sig0 ), length (n) ) / 
count  =0; 

for  ctrE  =  1 : length (E_sigy) 

for  ctrS  =  1 : length ( sigy_sig0 ) 
for  ctrN  =  1: length (n) 


nn  =n ( ctrN) ; 

sigy  =  sigy_sig0 (ctrS) *sig0 ; 

E  =  E_sigy (ctrE) *sigy_sig0 (ctrS) *sig0 ; 
options  =  optimset (' TolFun ', le-12 ,  '  Display ' ,  1  of f ' ) 
sigRl (ctrE, ctrS, ctrN)  =  fsolve(@(x)  ( sigy* (E/sigy* 
(x/E+  eRl)  ) Ann  -x) ,  sigO , options ) ; 

sigR2 (ctrE, ctrS, ctrN)  =  fsolve(@(x)  (sigy* (E/sigy* 
(x/E+  eR2)  ) Ann  -x) ,  sigO , options ) ; 

Ep  (ctrE, ctrS, ctrN)  =  ER_sigy(ctrE); 

Sp (ctrE, ctrS, ctrN)  =  sigy_sig0 (ctrS) ; 
np  (ctrE, ctrS, ctrN)  =  n(ctrN); 
count=  count+1 

end 

end 

end 

save  'matrix  redo2 ' 


else 

%  loads  the  resfence  stress  matrix  and  continues  the  program 
%  load  ' matrix_f ine4 ' 

%  load  ' matrix_coarse ' 

load  ' matrix_redo2 ' 

end 

checkl  =  sigRl(4,6,2)  /  (sigRl (4, 6, 2)  + 

E_sigy ( 4 ) *sigy_sig0 ( 6) *sigO*eRl ) An (2 ) . . . 

-  sigR2 ( 4 , 6,2)  /  ( sigR2 ( 4 , 6 , 2 )  + 

E_sigy ( 4 ) *sigy_sig0 (6) *sigO*eR2 ) An (2 ) 
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[Ep,Sp,np]  =  meshgrid (sigy_sigO , ER_sigy, n) ; 


o 

o 


ml  =  Ep . *Sp . *sigO . /sigRl ; 
m2  =  Ep . *Sp . *sigO . /sigR2 ; 

%Evaluate  the  fitting  functions 
Al  =  3.66556  +  0 . 024417 9*np; 

A2  =  6.06122-2 . 1 5  8  9 1 *  np ; 
q  =  29.0856  -  24.3547*np; 
p  =  1. 31861-0. 154675*np; 

kl  =  1.001  +0 . 2 610*np  -  0.5217*np.A2  +  0 . 1547*np . A3 ; 

k2  =  1.002  +0 . 7 637 *np  -  1.9200*np.A2  +  1.255*np.A3; 

jl  =  32.77  -  52 . 5 9*log (ml )  +  33 . 4 6* ( log (ml ) )  . A2  -4 . 8* ( log (ml ) )  . A3  .. 
+0.2147* (log (ml) ) . A4; 

j 2  -  8.817  -  12 . 7  3*log (ml )  +  11 . 99* ( log (ml ) )  . A2  -2 . 032 * ( log (ml ) )  . A3 
+0.1049* (log (ml) ) . A4; 

g  =  A2  +  (A1-A2 ) . / (1+ (m2 . /q) . Ap) ; 
fpl  =  kl . * j 1 ; 

fp2  =  k2  .  *  j  2  ; 

%  Calculate  Errors 

el  =  Cl. /sigRl  -  fpl; 

e2  =  C2./sigR2  -  fp2; 

e3  =  S . / (Ep . *Sp . *sig0 . *del2 )  -  g; 

e  =  abs(el)  +abs(e2)  +  abs(e3); 

%Find  location  of  values 

[r,c,u]  =  ind2sub ( size (e) , find (e  ==  min (min (min (e) ) )  )  ); 

%Output  new  values  of  E,  sig_y,  and  n 
new_E  =  Ep  (r , c, u) *Sp (r , c, u) *sig0* ( l-vA2 ) 
new_n  =  np (r, c,u) 
new_sigy  =  Sp  (r , c, u) *sig0 

min_el  =  min (abs (el ( : ) ) ) 
min_e2  =  min (abs  (e2  (  :  )  )  ) 
min_e3  =  min (abs  (e3  (  :  ) ) ) 
min_e  =  min (abs (e ( : ) ) ) 

Pd  =  abs (round ( [ (new_E-210e9) /210e9 ,  (new_n- .127)/ .127^  (new_sigy- 
400e6) /400e6] *100) ) 
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Spherical  (Hyperlastic) 


file  =  { 1 mu8sphere3 . dat ' ; ' mu8sphere3-reverse . dat ' } 


for  ctr_l  =  1 : length  ( file) 

fid=f open ( file { ctr_l },' r ') ;  %  Opens  .dat  file  for  Reading 

ctr  =  1; 
ctr2  =  1; 
ctr3  =  1; 
while  ~feof (fid) 

tline=fgetl (fid) ; 

if  isempty ( strf ind ( tline , ' N  ODE  OUTPUT'))  ==  0  %Scan 
till  header 

for  i=l : 9 

tline=f getl ( f id) ;  %  Read  10  lines  of  junk  lines 

end 

dataLine=f getl ( fid) ; %Grabs  the  reaction  force  and 

displacement 

data  =  sscanf (dataLine , ' %i  %f  %f'); 

F (ctr) =-data (2 ) /10A6; %  extracts  force  data  and  converts  to 
N 

d (ctr) =-data (3) /10A6; %  extracts  disp  data  and  converts  to  m 
ctr  =  ctr+1; 

end 

if  isempty ( strf ind (tline, 'C  ONTACT  OUTPUT'))  ==  0 
%Scan  till  header 

for  1*1:11 

tline=f getl ( fid) ;  %  Read  10  lines  of  junk  lines 

end 

dataLine=f getl ( fid) ; %Grabs  the  contact  stress 
if  dataLine  >0 

%  data  =  sscanf (dataLine,  ' %i  %s  %f  % f ' ) 
data  =  sscanf (dataLine, ' %f ') ; 

CA (ctr2) =data (1) / ( 1 0 A 6 ) A2;  %  extracts  contact  area  data 
and  converts  to  mA2 
end 

ctr2  =  ctr2+l; 

end 

if  isempty ( strf ind (tline, 'E  NERGY  OUTPUT'))  ==  0 
%Scan  till  header 

for  i=l : 22 

tline=f getl ( fid) ;  %  Read  10  lines  of  junk  lines 

end 

dataLine=f getl ( fid) ; %Grabs  the  contact  stress 
if  dataLine  >0 

%  data  =  sscanf (dataLine,  ' %i  %s  %f  % f ' ) 
data  =  sscanf (dataLine, ' %s  %s  %s  %s  %s  %f'); 

E (ctr3) =data (end) / ( 1 0 A 6 ) A2 ;  %  extracts  strain  energy 
data  and  converts  to  mA2 
end 
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ctr3  =  ctr3+l; 

end 

end 


matrix . F { ctr_l , : } 
matrix . d{ ctr_l , : } 
matrix . CA{ ctr_l , : 
matrix . E { ctr_l , : } 

end 

o  o 
o  o 


=  [0,  F (find (d>=0 ) ) ] 
=  [0,  d (find (d>=0 ) ) ] 
}  =  [0,  CA (find (d>=0 ) ) 
=  [0,  E (find (d>=0 ) ) ] 


%  Force  Displacement  plot  and  analysis 
R  =  1 2  0  e  -  6  ; 

h  =  linspace ( 0 , 0 . 4*R, 100 ) ; 

E  =  3*7 . 9  7  e  3 ; 
v  =  0.499; 

ER  =  E/ ( l-vA2 ) ; 


%N 

%m 

%mA2 

o 

o 


alpha  =  20; 

ac  =  sqrt (R*h) ; 

P  =  4/3*ER*RA.5*h. Al. 5; 


plot  (h,P,  '  — ') 
hold  on 

gr_str  =  [ 1 r+ 1 ; 1 kx 1 ;  ' ko ’ ; ’ bs ' ; ’ mo ' ] ; 

for  ctr  =  1 : length ( file) 

plot (matrix . d{ ctr ,  : } , matrix . F{ ctr ,  : } ,  gr_str (ctr, : ) ) 

end 

grid  on 

xlabel (' Displacement ,  h') 
ylabel ( ' Force,  N ' ) 
grid  on 

legend ( ’ Analytic ' , f ile { 1 } ,  ’ Location  1 ,  1  Northwest 1 ) 

%  legend ( 1 \alpha  =  97.1100  ’ , 1 \alpha  =  61.9200  1 , ' \alpha  =  59.2200 
’  ,  1  Location ' ,  ' Northwest ' ) 
figurehandle  =  gcf; 

set ( f indall (figurehandle,  ’ type  f ,  1  text ' ) ,  ' font size ' , 14 ,  ’ f ontweight ’ ,  ’ bol 
d’) 

%  Fitting  and  plotting  the  stress  strain  curve 


figure 
hold  on 

for  ctr  =  1 : length ( file) 

Fs  =  matrix . F { ctr } 

CAs  =  matrix . CA{ ctr } 
ds  =  matrix . d{ ctr } 

matrix . strain { ctr ,: }  =  [0,  0 . 2 *sqrt (ds ( 2 : end)  . /R)  ]; 

matrix . stress { ctr ,: }  =  [0,  Fs(2:end)  ./  CAs(2:end)] 

plot (matrix . strain { ctr , : } , matrix . stress { ctr , : } , gr_str (ctr, : ) ) 
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init  =  diff (matrix . stress { ctr, :}). /diff (matrix . strain { ctr, :} ) 
initE  =  init ( 1 ) 


s  =  fitoptions ( 'Method' , ' NonlinearLeastSquares ' ,  . . . 

' Robust ' , ' Bisquare '  ,  . . . 

' Lower ',[0,0],... 

'Upper' ,[],... 

' St art Point ' ,  [initE *40/3 /pi/0. 75,5] , . . . 

' Dif fMinChange ' , le-8 ,  .  .  . 

' Dif fMaxChange ' ,  . 1 ,  . . . 

'Maxlter '  ,  1000,  .  .  . 

' TolFun ' , eps , . . . 

' TolX ' ,  eps ) ; 

form  =  [  ' B/a* (  -(l-x).A(a-l  )+  (l-x).A  (-a/2-1))1]; 

f  =  f ittype ( f orm,  'independent',  'x',  'coefficients',  {'B','a'}, 

' options ' , s ) ; 

[fitobj,gof]  =  fit (matrix . strain { ctr ,:}' ,  matrix . stress { ctr ,:}' , 
f) 

cf (ctr,  :)  =  coef fvalues ( f itob j ) ; 

plot (matrix .strain {ctr,  : },  cf  ( ctr , 1 ) /cf (ctr, 2) * . .  . 

(  -  ( 1 -matrix . strain { ctr ,  : } )  . A  ( cf ( ctr , 2 ) -1 )  .  .  . 

+  (1-matrix. strain{ctr,  :}). A  (-cf (ctr, 2) / 2  —  1 ) ),' k'  ); 

end 

mu0_out_sig  =  cf ( : , 1 ) *3/40*pi* ( l-vA2 ) 
alpha_out_sig  =  cf(:,2) 
hold  on 

plot (matrix. strain { 1,  : } , matrix. strain { 1,  : } *ER*20/3/pi)  %plots  the 
%linear  elastic  solution 
grid  on 

xlabel ( ' Strain ' ) 
ylabel (' Stress ,  Pa') 
grid  on 

legend (' \mu_0  =  7.97kPa,  \alpha  =  20', 'Reverse  Analysis ',' Linear 
Elastic '  ,  ' Location ' ,  ' Northwest ' ) 

%  legend (' \alpha  =  97.1100  ' , 1 \alpha  =  61.9200  ' , ' \alpha  =  59.2200 
'  ,  ' Location ' ,  ' Northwest ' ) 
figurehandle  =  gcf; 

set ( f indall (figurehandle, ' type ' , ' text ' ) , ' font size ' , 14 , ' f ontweight ' , ' bol 
d '  ) 
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Appendix  C:  Hawkmoth  Rearing  [21] 


The  AFIT  Flapping  wing  MAV  research  group  receives  hawkmoth  pupae  on  a 
regular  basis  from  Dr.  Mark  Willis  at  Case  Western  University.  Dr.  Willis’  lab  contains  a 
thriving  colony  of  Manduca  Sexta,  which  produces  scores  of  moths  each  week.  The  most 
challenging  part  of  raising  the  hawkmoths,  hatching  the  eggs  and  feeding  the  caterpillars, 
is  already  complete  when  we  receive  the  pupae  in  the  mail.  All  they  need  are  a  proper 
light  cycle  and  the  right  temperature  in  order  to  eclose  (hatch  into  adulthood).  Figure  63 
shows  a  hawkmoth  pupa.  The  specimen  is  on  its  back  with  the  head  pointed  toward  the 
left.  The  right  forewing  can  be  seen  wrapped  around  midsection  of  the  body.  The 
abdomen  with  its  many  segments  and  spiracles,  points  to  the  right.  The  “handle”  object 
protruding  from  the  head  is  the  proboscis,  folded  several  times.  When  extended,  the 
proboscis  of  the  adult  M. sexta  can  be  as  long  as  it’s  body,  and  is  used  for  feeding  on 
flower  nectar  while  hovering. 


■L!*  proboscis  __ 

PI*-  win? 


Figure  63:  Typical  M.  sexta  pupa  and  diagram  of  individual  parts 
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A  lxlxl  foot  cubic  terrarium  holds  the  pupae  prior  to  eclosion.  The  bottom  of  the 
terrarium  is  lined  with  approximately  1  inch  of  wood  shavings  beneath  a  layer  of  paper 
towels.  This,  along  with  frequent  cleaning,  is  necessary  because  with  each  moth’s 
eclosion  comes  a  release  of  large  quantities  of  waste  ( meconium )  which  is  the  by-product 
of  the  transformation  from  caterpillar  to  moth  (Reinecke,  Buckner  and  Grugel  1980).  The 
front  of  the  terrarium  has  doors  which  can  swing  open  for  access  from  the  front  and  the 
entire  top  glass  pane  can  be  removed  for  access  from  the  top.  The  most  important  feature 
is  the  back  wall  made  of  textured  foam.  This  wall  gives  the  freshly  eclosed  moth  a 
surface  to  climb  up,  which  is  an  absolute  necessity.  The  moth  must  climb  off  of  the 
ground  in  order  to  pump  fluids  through  the  veins  in  its  wings  to  stretch  them  out  before 
they  harden.  Typically,  the  moth  finds  a  position  on  the  wall  in  about  10  minutes  and  has 
fully  inflated  its  wings  20  minutes  later.  Figure  64  shows  two  young  adult  hawkmoths 
which  hatched  only  a  few  minutes  apart. 


Figure  64:  Two  freshly-enclosed  M.  Sexta  specimens.  One  is  only  10  minutes  old  and 

has  not  inflated  its  wings. 
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M.sexta  thrives  in  the  warm  temperatures  of  the  southern  North  American 
summer.  The  pupae  prefer  a  summertime  light  cycle  of  long  days  and  short  nights,  as 
well  as  warm  summertime  temperatures.  The  light  cycle  for  the  AFIT  moths  has  been  set 
to  14  hours  of  light  and  10  hours  of  darkness  (Willis  2011).  The  light  cycle  is 
accomplished  using  a  standard  outlet  timer  wired  to  a  string  of  LED  lights  secured  around 
the  inner  walls  of  a  cardboard  box  (Figure  65).  The  terrarium  with  the  pupae  is  placed 
within  this  box  and  the  lid  is  then  closed,  allowing  no  light  from  the  outside.  The 
temperature  of  the  vivarium  is  set  to  80  degrees  Fahrenheit  and  the  humidity  is  kept  at  40 
percent  in  order  to  replicate  summertime  conditions. 


Figure  65:  Enclosure  for  the  pupae:  (Left)  Closed  box  creates  day/night  conditions 
with  a  tinier  controlling  the  LED  lights,  shown  in  part  (Right) 

The  terrarium  is  checked  daily  for  newly  emerged  adult  moths.  They  tend  to 
eclose  at  “dusk,”  or  shortly  after  their  light  cycle  switches  from  light  to  dark.  Since  the 
cardboard  box  enclosure  around  the  terrarium  allows  this  light  cycle  to  be  set  for  any 
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time  of  day,  the  light-to-dark  transition  has  been  set  to  occur  mid-afternoon.  That  way, 
the  moths  will  eclose  in  the  afternoon  and  the  terrarium  can  be  checked  for  adults  when 
heading  home  for  the  day.  When  an  adult  moth  is  found  in  the  pupa  terrarium  during 
daily  inspections,  it  is  transferred  to  a  mesh  cage  which  is  constantly  open  to  the  light  of 
the  vivarium.  M.sexta  is  a  nocturnal  species,  and  as  such  is  inactive  during  daylight.  The 
adults  are  therefore  docile  and  essentially  dormant  as  long  as  the  lights  are  kept  on.  If  the 
lights  are  switched  off,  however,  the  adults  will  fly  for  hours  against  the  walls  of  the 
enclosure  and  damage  their  wings.  This  is  to  be  avoided  because  much  of  the  research 
that  goes  on  among  the  AFIT  FWMAY  research  group  requires  intact  forewings. 
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